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OBJECTIVE 


To  improve  MINIMUF-3.5. 


RESULTS 


1.  The  model  was  improved  by  comparing  the  predicted  data  to  fQF2 
values  measured  at  30  sites. 

2.  The  sunspot  number  variation  in  MINIMUF-3.5  was  deleted  and  the 
constant  58.0  called  A1  in  MINIMUF-3.5  was  replaced  by  a  linear  equation  in 
sunspot  number. 

3.  The  M- factor  portion  of  the  algorithm  was  modified  to  include 
sunspot  number,  seasonal,  and  diurnal  variations  using  over  7200  observed 
oblique  sounder  median  MOFs  on  39  paths  in  the  fitting  process. 

A.  A  special  fQF2  model  was  added  for  use  in  the  polar  regions. 

RECGM1ENDATI0NS 

1.  Use  MINIMUF-85  in  all  applications  now  using  MINIMUF-3.5. 

2.  Use  the  fQF2  model  in  MINIMUF-85  in  ray-tracing  applications  where 

f  F2  from  MINIMUF-3.5  is  now  used, 
o 

3.  Use  the  M-factor  model  in  MINIMUF-85  to  determine  a  M3000  factor  for 
obtaining  mirror  reflection  heights  for  oblique  incidence  propagation  rather 
than  a  fixed  F2-region  height. 

A.  Use  the  resultant  value  of  mirror  reflection  height  to  determine 
take-off  angles  for  use  in  determining  antenna  gains  and  path  loss  in 
applications  where  now  a  fixed  value  of  take-off  angle  at  a  given  range  is 
used. 


1 


5.  Use  the  f  F2  data  base  to  improve  the  f  F2  representation  of 
o  ° 

MINIMUF-85  by  determining  the  sunspot  number,  seasonal,  and  geographic 
dependencies  of  the  parameters  Aq,  and  cos 

6.  Introduce  the  effects  of  underlying  layers  on  M-factor  estimation. 

7.  Continue  to  improve  the  polar  representation. 
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INTRODUCTION 


The  effective  operation  of  long  distance  high  frequency  (HF)  systems  is 
increased  in  proportion  to  its  ability  to  predict  variations  in  the  iono¬ 
sphere,  since  such  ability  permits  the  selection  of  optimum  frequencies, 
antennas,  and  other  circuit  parameters.  Most  variations  in  HF  system 
performance  are  directly  related  to  changes  in  the  ionosphere,  which  in  turn 
are  affected  in  a  complex  manner  by  solar  activity,  seasonal  and  diurnal 
variations  as  well  as  latitude  and  longitude. 

Originally,  manual  methods  were  developed  for  analyzing  these  effects  on 
HF  circuits  of  short,  intermediate, and  long  distances  (reference  1).  Because 
the  manual  methods  are  laborious  and  time  consuming,  various  organizations 
developed  computer  programs  to  analyze  HF  circuit  performance.  A  commonly 
predicted  parameter  in  these  programs  is  the  maximum  usable  frequency  (MUF). 
The  MUF  is  the  highest  frequency  that  can  be  propagated  by  ionospheric 
refraction  between  given  points  at  a  given  time. 

However,  these  computer  programs  depended  on  the  use  of  main  frame 
computers.  In  the  day-to-day  management  of  frequency  assets  of  a  communica¬ 
tion  system,  the  usual  procedure  was  to  rely  on  a  set  of  so-called  "propaga¬ 
tion  charts"  produced  by  these  programs.  To  meet  publication  deadlines,  these 
charts  had  to  be  produced  months  in  advance.  Hence,  any  possibility  of  real¬ 
time  prediction  was  not  possible. 

In  1978  a  simple  semiempirical  algorithm,  called  Mj.NIMUF-3.5,  was 
developed  for  predicting  the  MUF  on  small  mobile  propagation  forecast 
(PROPHET)  terminals  (references  2,  3).  With  this  tool  it  was  possible  to 
develop  a  variety  of  new  forecast  applications  where  the  use  of  the  large- 
scale  propagation  programs  in  the  operational  environment  was  not  practical. 

In  its  initial  development,  verification  of  MINIMUF-3.5  was  done  by 
comparing  the  predictions  with  oblique- incidence  sounder  data.  The  data  base 
used  encompassed  196  path  months  (4704  test  points)  of  observed  maximum 
observed  frequencies  (MOF)  over  23  different  HF  sounder  paths.  MINIMUF-3.5 
was  found  to  have  an  rms  error  of  3.8  MHz. 
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Subsequent  to  this  comparison,  a  more  complete  verification  of  MINIMUF 
3.5  was  made  (reference  A).  In  this  test,  4668  MOFs  measured  on  25  paths  are 
compared  against  the  predicted  values  from  ITSA-1  (reference  5)  and  HFMUFES  4 
(references  6,  7).  The  data  were  divided  into  subsets  to  see  the  effect  of 
particular  paths,  path  length  and  orientation,  season,  month,  sunspot  number, 
diurnal  trends,  geographic  region  and  sounder  type.  MINIMUF-3.5  had  a  bias  of 
0.08  MHz  low  (0.6  percent  low)  and  an  rms  error  of  3.7  MHz  (3.6  percent), 
was  least  accurate  during  the  sunrise  and  sunset  transition  hours  and  for  path 
lengths  5000  to  7000  km  than  it  was  for  other  times  of  the  day  and  path 
lengths.  Linear  regression  analysis  showed  that  its  errors  in  predicted  MUF 
at  sunrise  and  for  path  length  lengths  5000  to  7000  km  are  nonlinear  and  could 
probably  be  attributed  to  the  M  factor  part  of  the  calculation. 


During  July  and  August  1982,  a  field  test  was  conducted  from  Ft.  Lewis, 
WA;  Ft.  Leavenworth,  KS;  and  Ft.  Knox,  KY;  to  Ft.  Bragg,  NC,  using  oblique- 
incidence  sounders.  One  of  the  models  used  in  the  test  for  operational 
frequency  selection  included  MINIMUF-3.5.  A  comparison  was  made  between  the 
predicted  MUFs  for  those  paths  and  the  observed  MOFs.  When  the  F-region  MOFs 
are  compared  against  the  predicted  MUFs  from  MINIMUF-3.5,  it  was  found  that 
MINIMUF-3.5  consistently  predicted  high.  This  experiment  was  conducted  during 
a  period  in  the  solar  cycle  for  which  solar  activity  was  unusually  high  and 
variable.  The  mean  sunspot  number  was  113,  its  standard  deviation  was  £8,  and 

the  peak  value  was  272. 


This  report  describes  the  development  of  an  improved  version  of  MINIMUF- 
3.5,  called  MINIMUF-85.  The  version  developed  to  predict  accurate  maximum 
useable  frequencies  (MUFs)  under  conditions  of  anomalously  high  sunspot 
numbers,  to  predict  foF2  values  suitable  for  ray-tracing  applications,  to 
predict  M  factor  values  useable  for  determining  the  mirror  height  of  reflec 
tion  for  oblique  incidence  propagation,  and  to  improve  its  accuracy  tor  paths 
into  or  crossing  the  polar  regions.  This  version  includes  sunspot  number 
dependence  in  both  the  f0F2  ana  the  M  factor  calculations.  The  M  factor 
portion  of  the  algorithm  also  was  modified  to  include  seasonal  and  diurnal 
variations.  Appendices  A  and  B  give  listings  of  MINIMUF-85  in  the  BASIC  and 

FORTRAN  languages,  respectively. 
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BACKGROUND 


MINIMUF  3.5  is  a  semi-empirical  model  developed  in  1978  (the  initial 
algorithm  was  called  MINIMUF-3)  to  provide  a  maximum -usable-  frequency  (MUF) 
prediction  capability  suitable  for  use  on  small  (micro)  computers  where  time 
and  storage  limitations  exist.  The  theory  and  method  used  in  the  development 
of  the  MINIMUF  3.5  algorithm  has  been  documented  in  several  earlier  reports 
and  will  not  be  presented  here  (references  2  and  3). 

The  expression  for  the  MUF  used  in  a  MINIMUF  3.5  is  given  by 

MUF  *  M  •  f  F2  (1) 

o 

where  M  is  the  obliquity,  or  M-,  factor  which  reflects  the  dependence  of  the 
MUF  on  transmission  path  length.  The  parameter  fQF2  is  the  critical  or 
penetration  frequency  at  vertical  incidence  for  the  F2  layer. 

In  particular,  we  have 

M  -  {1  +  2.5[sin(2.54'?)]3/2}  •  Gj  •  Gj  •  G3  (2) 

where  is  the  minimum  great  circle  distance  between  transmitter  and  receiver. 
The  various  constants  in  the  bracketed  term  in  equation  2  are  determined  by 
fitting  this  expression,  without  the  G^  i  *1,2, 3, to  an  exact  transmission 
curve  for  a  parabolic  layer  height  of  290  km  and  a  ratio  of  height  of  maximum 
of  electron  density  to  half -width  of  the  F2  layer  of  0.4  (reference  2).  The 
multipliers  G^.  provide  small  corrections  to  the  MUF  for  known  systematic 
departures  from  the  median  behavior  under  certain  conditions  of  path  geography 
or  season  (reference  2). 

The  expression  for  the  critical  frequency  used  in  MINIMUF  3.5  is 

V2-(‘  +  r)[*o  +  *i^^l/2  <3) 

where  Rq  Aq,  are  constants  and  R  is  the  12  month  running  mean  sunspot 
number.  The  constants  in  equation  3  were  determined  by  iteratively  adjusting 
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the  model  in  a  "real-time"  mode,  to  36  path  months  of  data  chosen  to  represent 
a  range  of  transmission  path  types  (reference  2). 

In  equation  3,  xeff  is  an  "effective"  solar  zenith  angle.  CosXeff  is 
modeled  as  the  logged  response  of  a  dynamic  linear  system,  "driven"  by  the 
instantaneous  value  of  cos  x*  By  using  an  effective  value  of  the  zenith 
angle,  recognition  is  given  to  the  fact  that  the  F2  layer,  unlike  the  E  and  D 
layers,  does  not  show  a  relatively  simple  cos  \  diurnal  dependence  on  x*  The 
dynamical  behavior  of  the  F2  layer  is  more  complicated  because  various  other 
dependencies  make  simple,  accurate  modeling  difficult.  In  kedping  with  the 
simplistic  nature  of  the  model,  defining  an  effective  x  allows  relatively 
accurate  modeling  without  explicitly  including  these  other  dependencies 
(reference  2). 
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DEVELOPMENT  PROCEDURE 


This  section  will  present  how  MINIMUF-85  is  developed  from  MINIMUF-3.5 
using  both  vertical  incidence  and  oblique  incidence  sounder  data.  Measure¬ 
ments  of  f  F2  data  at  30  selected  sites  were  used  to  improve  the  f  F2  model, 
o  o 

Measurements  on  39  paths  were  collected  to  form  a  data  base  of  over  7200 
observed  oblique  sounder  median  MOFs.  These  data  were  used  to  determine  the 
optimum  location  of  control  points  to  be  used  to  introduce  a  geomagnetic 
latitude  dependence,  to  improve  the  M-factor  model,  and  to  verify  the  polar 
region  model.  To  compare  the  model  against  the  measured  data  sets  a  data 
screening  program  is  used. 

VERTICAL  INCIDENCE  SOUNDER  DATA  BASE  (f  F2) 

o 

The  derivation  of  the  critical  frequency  expression  used  in  the  new  model 
to  be  described  in  this  report  required  the  construction  of  an  fQF2  data  base. 
For  this  purpose  we  were  fortunate  to  have  access  to  the  complete  collection 
of  hourly  monthly  media  f  F2  data,  collected  since  1930  at  many  vertical 
incidence  sounder  sites  and  reported  to  the  World  Data  Center  in  Boulder, 
CO.  With  access  to  this  large  amount  of  data  we  were  able  to  be  selective  in 
determining  which  sites  and  time  periods  to  include  in  on1'  data  base.  The 
selection  criteria  used  are  based  on  geographic  distribution  of  the  sites  and 
availability  of  data  covering  the  range  of  sunspot  numbers  we  wish  to  include. 

Figure  1  shows  the  geographic  distribution  of  the  30  selected  sites. 
Table  1  gives  the  name  and  location  included  for  each  site  in  the  fQF2 
data  base.  With  the  available  data  we  attempted  to  obtain  maximum  distribu¬ 
tion  over  all  the  continents  in  order  to  include  all  known  geographical 
effects  on  the  fQF2  in  the  data.  The  fact  that  most  of  the  longer  operating 
stations  reporting  data  are  in  mid- latitudes  (-20° to 60°)  means  that  this  region 
of  the  world  is  usually  over-represented  in  comparison  with  other  regions, 
i.e.,  polar  and  equatorial  regions.  This  is  perhaps  the  case  in  our  fQF2  data 
base  since  the  second  selection  criterion ,  data  covering  a  wide  range  of  sun¬ 
spot  numbers,  forced  the  inclusion  of  many  of  these  sites. 
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30-  WEST  0"  EAST  30' 


ISO' EAST  180’ WEST  ISO'  120‘  SO' 


Figure  1 .  Geographic  distribution  of  sites  in  f0F2  data  base. 


Table  1.  Name  and  geographic  location  of  sites  in  fQF2  data  base. 


SITE 

Latitude 

Longitude 

SITE 

Latitude 

Longitude 

Resolute  Bay 

74. 7N 

265. IE 

Washington,  DC 

37. 4N 

237. 8E 

Murmansk 

69.0 

33.0 

San  Francisco 

37.4 

237.8 

College 

64.9 

212.2 

White  Sands 

32.3 

253.5 

Anchorage 

61.2 

210.1 

Ahmedabad 

23. 

72.6 

Kjeller 

60.0 

11.1 

Maui 

20.8 

203.5 

Churchill 

58.8 

265.8 

Mexico  City 

19.4 

260.3 

Moscow 

55.5 

37.3 

Puerto  Rico 

18.5 

292.8 

Juliusruh/Rugen 

54.6 

13.4 

Ibadan 

7.4 

3.9 

Lindau 

51.6 

10.1 

Paramaribo 

5.8 

304.8 

Slough 

51.5 

359.4 

Singapore 

1.3 

103.8 

Dourbes 

50.1 

4.6 

Kinshasa-Binza 

4.5S 

15.2 

Winnipeg 

49.8 

265.6 

Huncayo 

12. OS 

284.7 

Freiburg 

48.1 

7.6 

Hobart 

42. 9S 

147.3 

St.  Johns 

47-6 

307.3 

Godley  Head 

43. 6S 

172.8 

Wakkanai 

45.4 

141.7 

Kerguellen 

49. 4S 

70.3 
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For  each  site  in  table  1,  the  data  base  contains  a  winter  month,  a  summer 
month,  and  a  spring  and  fall  equinox  month  in  each  of  four  sunspot  number 
ranges,  low  (0  -  30),  medium  (31  -  100),  high  (101  -  150)  and  very  high 

(>150).  This  gives  a  total  of  16  months,  each  with  24  hourly  median  data 

values,  at  each  site.  In  total  then,  there  are  11,520  data  points  in  the 
data  base.  The  months  and  years  for  each  site  are  given  in  table  2. 

We  feel  that  this  represents  a  reasonably  complete  and  inclusive 

description  of  the  geographical  and  solar  dependencies  of  the  fQF2  parameter. 

OBLIQUE  INCIDENCE  SOUNDER  DATA  BASE  (HOF) 

The  oblique  sounder  data  base  that  was  assembled  was  derived  from  a 
variety  of  sources  and  spans  the  period  between  1958  and  1976.  This 

data  base  represents  over  one  complete  solar  sunspot  cycle  of  propagation 
data.  Attempts  were  made  to  make  the  data  base  as  diverse  as  possible, 
including  a  variety  of  different  path  lengths,  orientations,  geographical 
locations,  and  sunspot  numbers. 

The  oblique  sounder  data  set  consists  of  304  path-months  of  median  hourly 
MOF  values  derived  from  39  different  HF  transmission  paths.  In  addition  to 
the  original  oblique  sounder  data  base  reported  by  Sailors  et  al.  (reference 
4),  this  data  base  included  14  additional  paths.  Table  3  summarizes  the  addi¬ 
tional  HF  oblique  sounder  data  base.  The  locations  of  all  the  paths  are  shown 
in  figure  2  except  for  paths  for  which  the  scale  is  too  small  to  illustrate. 

Of  the  39  paths,  the  longest  path  was  7808  km  and  the  shortest  path  was  192 

km.  Table  4  shows  the  percentage  of  the  sample  in  each  path  length  range. 
Table  5  shows  the  percentage  of  the  sample  in  different  path  orientation 
categories.  The  north-south  paths  are  those  which  lie  nominally  within  il5° 

of  0°  or  180°  bearing.  The  east-west  paths  are  those  which  fall  nominally 

within  ±15°  of  a  90°  or  273°  bearing.  The  paths  which  did  not  meet  either 
criterion  were  put  in  the  "other"  category.  Table  6  shows  the  percentage  of 
the  sample  categorized  according  to  the  geomagnetic  latitude  location  of 
control  points.  Table  7  shows  the  percentage  of  the  sample  in  each  sunspot 
number  (SSN)  category. 
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Table  2.  f  F2  data  at  each  site, 
o 


Site 

Summer 

Fall 

Winter 

Spring 

6/54 

6/60 

9/54 

9/70 

12/53 

12/59 

4/54 

4/60 

Resolute  Bay 

6/66 

6/57 

4/62 

4/58 

12/61 

12/57 

4/62 

4/58 

6/64 

6/68 

9/64 

9/68 

12/63 

12/68 

4/65 

4/69 

Murmansk 

6/66 

6/59 

9/66 

9/58 

12/61 

12/58 

4/62 

4/59 

6/54 

6/60 

9/54 

9/46 

12/53 

12/59 

4/54 

4/60 

College 

6/52 

6/57 

9/50 

9/58 

12/50 

12/57 

4/42 

4/58 

6/54 

6/60 

9/54 

9/60 

12/53 

12/59 

4/54 

4/60 

Anchorage 

6/62 

6/57 

9/62 

9/58 

12/61 

12/57 

4/62 

4/58 

6/54 

6/56 

9/54 

9/49 

12/54 

12/49 

4/54 

4/53 

Kjeller 

6/55 

6/57 

9/55 

9/57 

12/55 

12/57 

4/50 

4/58 

6/54 

6/60 

9/54 

9/70 

12/53 

12/59 

4/54 

4/60 

Churchill 

6/66 

6/57 

9/62 

9/58 

12/61 

12/57 

4/62 

4/58 

6/54 

6/60 

9/54 

9/70 

12/53 

12/59 

4/54 

4/60 

Moscow 

6/66 

6/57 

9/62 

9/58 

12/61 

12/57 

4/62 

4/58 

6/64 

6/68 

9/64 

9/68 

12/63 

12/68 

4/65 

4/69 

Juliusruh/Rugen 

6/66 

6/59 

9/66 

9/58 

12/61 

12/58 

4/62 

4/59 

6/54 

6/60 

9/54 

9/68 

12/53 

12/59 

4/54 

4/60 

Lindau 

6/66 

6/57 

9/62 

9/58 

12/61 

12/57 

4/62 

4/58 

6/54 

6/60 

9/54 

9/70 

12/53 

12/59 

4/54 

4/60 

Slough 

6/66 

6/57 

9/62 

9/58 

12/61 

12/57 

4/62 

4/58 

6/64 

6/68 

9/64 

9/68 

12/63 

12/68 

4/65 

4/69 

Dourbes 

6/66 

6/59 

9/66 

9/58 

12/61 

12/58 

4/62 

4/59 

6/54 

6/60 

9/54 

9/70 

12/53 

12/59 

4/54 

4/60 

Winnipeg 

6/66 

6/57 

9/62 

9/58 

12/72 

12/57 

4/62 

4/58 

6/54 

6/60 

9/54 

9/70 

12/53 

12/59 

4/54 

4/60 

Freiburg 

6/66 

6/57 

9/62 

9/58 

12/61 

12/57 

4/62 

4/58 

6/54 

6/60 

9/64 

9/70 

12/53 

12/59 

4/54 

4/60 

St.  Johns 

6/66 

6/57 

9/62 

9/58 

12/61 

12/57 

4/62 

4/58 

6/54 

6/60 

9/54 

9/70 

12/54 

12/59 

4/54 

4/60 

Wakkanhi 

6/66 

6/57 

9/62 

9/58 

12/61 

12/57 

4/62 

4/58 

6/54 

6/60 

9/54 

9/60 

12/53 

12/59 

4/54 

4/60 

Washington,  DC 

6/66 

6/57 

9/62 

9/58 

12/61 

12/57 

4/62 

4/58 
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Table  2  (continued). 


Site 

Sumner 

Fall 

Winter 

6/52 

6/57 

9/55 

9/57 

12/50 

12/578 

San  Francisco 

6/52 

6/57 

9/55 

9/57 

12/50 

12/57 

6/54 

6/60 

9/54 

9/70 

12/53 

12/59 

White  Sands 

6/66 

6/57 

9/62 

9/58 

12/61 

12/57 

6/64 

6/68 

9/64 

9/68 

12/63 

12/68 

Ahmedabad 

6/66 

6/59 

9/66 

9/58 

12/61 

12/58 

6/54 

6/60 

9/54 

9/70 

12/53 

12/59 

Maui 

6/66 

6/57 

9/62 

9/58 

12/61 

12/57 

6/64 

6/68 

9/64 

9/68 

12/63 

12/68 

Mexico  City 

6/66 

6/59 

9/66 

9/59 

12/61 

12/58 

6/44 

6/51 

9/54 

9/46 

12/53 

12/49 

Puerto  Rico 

6/45 

6/57 

9/50 

9/57 

12/51 

12/57 

6/54 

6/60 

9/54 

9/70 

12/54 

12/59 

Ibadan 

6/66 

6/57 

9/62 

9/53 

12/61 

12/57 

6/64 

6/60 

9/64 

9/60 

12/64 

12/60 

Paramaribo 

6/62 

6/58 

9/61 

9/52 

12/61 

12/57 

6/54 

6/60 

9/54 

9/60 

12/53 

12/59 

Singapore 

6/62 

6/57 

9/62 

9/58 

12/61 

12/57 

6/543 

6/60 

9/54 

9/60 

12/53 

12/59 

Kinshasa-Binza 

6/55 

6/57 

9/55 

9/58 

12/61 

12/57 

6/54 

6/60 

9/54 

9/70 

12/53 

12/59 

Huncayo 

6/66 

6/57 

9/62 

9/58 

12/61 

12/57 

6/54 

6/51 

9/54 

9/46 

12/53 

12/47 

Hobart 

6/52 

6/57 

9/50 

9/58 

12/50 

12/57 

6/54 

6/60 

9/54 

9/70 

12/53 

12/59 

Godley  Head 

6/66 

6/57 

9/62 

9/58 

12/61 

12/57 

6/64 

6/68 

9/64 

9/68 

12/63 

12/68 

Kerguellen 

6/66 

6/59 

9/66 

9/58 

12/61 

12/58 

Spring 


4/45  4/48 
4/45  4/48 

4/54  4/60 

4/62  4/58 

4/65  4/69 

4/62  4/59 

4/54  4/60 

4/62  4/58 

4/65  4/69 
4/62  4/59 

4/44  4/50 

4/45  4/48 

4/54  4/60 

4/62  4/58 

4/65  4/60 

4/62  4/58 

4/54  4/60 

4/62  4/58 

4/54  4/60 

4/61  4/58 

4/54  4/60 

4/62  4/58 

4/54  4/50 

4/52  4/58 

4/54  4/60 

4/62  4/58 

4/65  4/69 
4/62  4/58 
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Table  3.  Additional  HF  propagation  oblique  sounder  data 
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Table  4.  Percentage  of  sample  oblique  data  in  each  path  length  range. 


Length 

Number  of  Hours 

Percentage  of  Sample 

L  <  1000 

216 

3.0 

1000  <  L  <  2000 

96 

1.3 

2000  <  L  <  3000 

1000 

13.7 

3000  <  L  <  4000 

1057 

14.5 

4000  <  L  <  5000 

1608 

22.1 

5000  <  L  <  6000 

1567 

21.5 

6000  <  L  <  7000 

1110 

15.3 

7000  <  L  <  8000 

622 

8.5 

Total 

7276 

i 

100.0 

Table  5.  Percentage  of  sample  oblique  data  in  path  orientation  categories. 


Path  Orientation 

Number  of  Hours 

Percentage  of  Sample 

North/South 

2629 

36.1 

East/West 

1107 

15.2 

Other 

3516 

48.3 

Table  6.  Percentage  of  sample  oblique  data  in  geomagnetic  latitude 
categories. 


Path  Type 

Number  of  Hours 

Percentage  of  Sample 

Transequatorial 

2031 

27.9 

Low  Latitude 

984 

13.5 

Mid-latitude 

2712 

37.3 

High  Latitude 

977 

13.4 

Transauroral 

572 

7.9 

Table  7.  Percentage  of  sample  in  each  sunspot  number  category. 


Sunspot  Number 
(Cycle  Phase) 

10-30  (minimum) 

Number  of  Hours 

Percentage  of  Sample 

1865 

25.6 

31-60  (rise  and  decline) 

1121 

15.4 

61-90  (near  maximum) 

1860 

25.6 

91-120  (maximum) 

2101 

28.9 

121-150  (high  maximum) 

257 

3.5 
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DATA  SCREENING 


In  the  comparison  of  a  program  against  data,  it  is  desirable  to  subdivide 
the  data  base  into  subsets  according  to  variables  influencing  the  predicted 
and  observed  results  (e.g.,  path  length,  season,  month,  geomagnetic  latitude, 
sunspot  number,  local  time  at  path  midpoint,  etc.)*  To  accomplish  this,  a 
computer  program  called  DASCR3  (acronym  for  data  screening  3)  is  used.  Each 
of  the  prediction  programs  is  run  for  each  of  the  paths  or  sites  in  the  data 
bases.  The  results  along  with  auxiliary  information  about  the  propagation 
situation  (e.g.,  path  length,  local  time  of  day,  sunspot  number,  etc.)  are 
stored  in  a  data  file  to  be  used  later  by  DASCR3. 

DASCR3 

DASCR3  is  a  program  designed  to  perform  data  screening  and  statistical 
comparison  on  two  large  matrices  of  observations.  For  each  set  of  matrices, 
up  to  10  sets  of  information  are  read  in  on  propositions  to  be  satisfied  and 
limits  on  a  selected  variable.  A  portion  of  each  matrix  is  read  in  and  tested 
for  each  set  of  propositions  in  turn.  For  each  subset  satisfying  a  given  set 
of  conditions,  the  variable  to  be  analyzed  is  stored  temporarily  on  disc.  The 
next  portion  of  each  matrix  is  then  read  in  and  screened  and  the  good  observa¬ 
tions  are  added  to  those  already  on  disc.  When  the  entire  matrix  has  been 
screened,  the  screened  data  are  then  read  into  core  and  the  difference  (or 
residual)  between  the  two  matrices  is  taken.  These  arrays  are  then  sorted  to 
ensure  maximum  computer  efficiency  for  the  statistical  evaluation.  Finally,  a 
statistical  evaluation  is  then  performed  of  the  screened  data  and  their 
residuals. 

An  example  of  the  output  from  DASCR3  is  given  in  figure  3.  In  this 
example,  MINIMUF  3.5(1)  is  compared  to  the  observed  data.  The  proposition  to 
be  satisfied  is  all  circuits  in  the  data  base  are  to  be  used.  The  variables 
being  compared  are  the  observed  MOF  and  predicted  MUF.  In  the  printout  the 
observed  data  are  represented  by  column  A  and  the  predicted  values  are  repre¬ 
sented  by  column  B.  The  residual  (the  observed  data  minus  the  predicted 
value)  is  given  by  column  D.  The  relative  residual  is  given  by  column  D/A, 
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Figure  3.  Example  of  output  from  DASCR3 


and  the  absolute  relative  residual  by  column  ABS(D)/A.  The  left  hand  side  of 
the  page  shows  the  statistics  calculated  for  each  of  these  columns.  In  addi¬ 
tion,  the  correlation  coefficient  between  the  observed  and  predicted  data  are 
given.  Included  also  are  the  slope,  intercept  and  mean  square  error  of  linear 
regression.  The  standard  error  of  the  estimate  of  Y  =  AX,  where  Y^  is  the 
measured  data  point  and  Xi  the  predicted  point,  is  given  as  well  as  the  slope 
A.  In  this  example,  7276  data  points  were  selected  by  DASCR3  from  7276  data 
points.  Note  that  the  average  absolute  relative  residual  for  this  case  is 
20.0  percent. 

SCREENING  DATA  BASE 

Each  version  of  the  computer  program  being  tested  is  run  to  produce  a 
data  base  corresponding  to  the  observed  data  base.  Auxiliary  information  out¬ 
putted  to  be  screened  included  universal  time  of  propagation,  month,  year, 
sunspot  number,  path  length  in  kilometers,  geographic  latitude  and  longitude 
of  the  path  midpoint,  the  local  time  at  the  path  midpoint,  the  path  orienta¬ 
tion  with  respect  to  north,  the  geomagnetic  latitude  at  each  of  the  control 
points,  the  predicted  MUF,  path  identification  number,  and  sounder  type. 

Before  the  actual  data  screening  has  begun,  data  points  in  both  observed 
and  predicted  bases  corresponding  to  observed  values  at  the  extremes  of  the 
particular  measuring  sounder  are  removed  from  the  data  base.  The  final  number 
of  hourly  values  in  the  data  base  was  7276  points. 

ANALYSIS  OF  RESIDUALS  BETWEEN  PREDICTIONS 
AND  OBSERVED  DATA 

An  indication  of  the  accuracy  of  the  numerical  predictions  of  MUF  can  be 
obtained  from  a  study  of  the  residuals  between  observed  data  and  predicted 
values.  The  terms  residual,  relative  residual,  and  absolute  relative  residual 
are  used  with  the  following  standard  meaning: 

residual  =  (observed  datum)  -  (predicted  value)  (A) 
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relative  residual 


residual _ 

observed  datum 


(5) 


absolute  residual 

absolute  relative  residual  -  observed  datum 


(6) 


Certain  statistical  measures  of  these  terms  have  proved  useful  in  past  iono¬ 
spheric  studies  in  comparing  predicted  and  observed  data  (reference  4).  These 
include 


(1)  The  average  residual  (av.  res.) 

(2)  Root  mean  square  residual  (rms  res.) 

(3)  The  mean  absolute  error  of  the  residual  (mae  res.) 

(4)  The  average  relative  residual  (av.  rel.  res.) 

(5)  The  root  mean  square  relative  residual  (rms  rel.  res.) 

(6)  The  mean  absolute  error  of  the  relative  residual  (mae  rel.  res.) 

(7)  The  average  absolute  relative  residual  (ave.  abs.  rel.  res.) 

(8)  Correlation  coefficient  between  observed  and  predicted  values 

(9)  The  standard  error  of  the  estimate  of  linear  regression 


Values  of  each  of  these  parameters  are  produced  by  DASCR3  as  can  be  seen  by 
examining  figure  3.  The  average  residual  and  the  average  relative  residual 
locate  the  center  of  the  distributions  of  error  and  are  sometimes  referred  to 
as  the  bias  in  the  estimate.  The  mean  absolute  errors  of  the  residual  and 
relative  residual  are  a  measure  range  of  the  error  and  are  the  first  moments 
about  the  average  residual  and  average  relative  residual,  respectively.  They 
provide  information  about  the  range  of  variation.  The  average  absolute  rela¬ 
tive  residual  is  a  measure  of  the  average  magnitude  of  the  error. 


The  root  mean  square  residual  and  relative  residuals  are  measures  of  the 
dispersion  in  the  error.  In  fact,  the  rms  residual  and  rms  relative  residual 
are  the  standard  deviations  of  the  error  about  the  origin  (zero  bias)  and  are 
related  to  the  standard  deviation  about  the  mean  according  to 

c2  -  v2  -  v,2  <»> 

where  v2  is  the  mean  square  error  (the  square  of  the  rms  error)  and  Vj  is  the 
bias.  When  the  bias  is  small  or  nearly  zero,  then  the  standard  deviation  and 
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the  rms  error  are  nearly  the  same.  Otherwise,  the  rms  error  is  larger  than 
the  standard  deviation. 

A  measure  of  the  degree  of  association  or  the  closeness  of  fit  between 
variables  is  given  by  the  correlation  coefficient.  The  degree  indicates  the 
strength  of  the  tendency  for  high  (or  low)  values  of  one  variable  to  be 
associated  with  high  (or  low)  values  of  the  other  variable. 


A  description  of  the  nature  of  the  relationship  between  variables  is 
called  regression  analysis  (reference  13).  Regression  analysis  is  concerned 
with  the  problem  of  describing  or  estimating  the  value  of  one  variable,  called 
the  dependent  variable,  on  the  basis  of  one  or  more  other  variables,  called 
independent  variables.  In  other  cases  regression  may  be  used  merely  to 
describe  the  relationship  between  known  values  of  two  or  more  variables. 


Regression  analysis  that  involves  the  determination  of  a  linear  relation¬ 
ship  between  two  variables  is  referred  to  as  simple  linear  regression.  Here, 
the  variable  y  is  given  as  y  =  a  +  bx  where  x  is  the  independent  variable  and 
y  is  the  dependent  variable.  The  coefficients  a  and  b  are  determined  in  the 
regression  analysis.  A  measure  of  the  success  of  linear  regression  analysis 
is  the  standard  error  of  the  estimate  given  by 

S 

y.x 

where  o  is  the  standard  deviation  in  the  observed  datum  and  y  is  the  corre- 

y 

lation  coefficient  between  the  observed  data  in  predicted  values.  If  the 
relationship  is  truly  linear,  then  the  bias  of  the  estimate  should  be  removed 
(or  made  nearly  zero).  An  estimate  of  the  standard  error  of  mean  is  given  by 


Sy-X  =  /n 


(10) 


A  measure  of  the  error  in  the  regression  coefficient  is  given  by 
1 1/2 


Sb  = 


JL2 i 


(11) 
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SUNSPOT  REPRESENTATION 


In  1984,  an  improved  version  of  MINIMUF-3.5  was  presented,  which  improved 
the  accuracy  of  the  algorithm  for  sunspot  numbers  greater  than  100  by  limiting 
the  growth  of  the  f  F2  portion  of  the  calculation  as  the  sunspot  number  became 
high  (reference  14).  It  retained  the  simplicity  of  MINIMUF-3.5  and  is  as 
accurate  as  MINIMUF-3.5  at  sunspot  numbers  less  than  100. 

For  years  of  high  solar  activity,  the  literature  indicates  that  an 
algorithm  with  a  limitation  on  the  increase  of  MUF  at  high  solar  activity  is 
desirable  (references  15,  16,  and  17).  Vasil'yeva  and  Kerblay  indicate  that 
there  are  three  types  of  dependencies  of  fQF2,  a  factor  in  the  MUF  computa 
tion,  on  solar  activity.  Type  I  includes  the  kinds  of  dependencies  that  can 
be  considered  as  linear  over  the  entire  range  of  sunspot  number  R.  Type  II 
behavior  is  characterized  by  fQF2  increasing  with  R  to  100-140  and  remaining 
constant  or  changing  insignificantly  with  further  increase  in  R.  In  type  III 
behavior,  the  dependence  is  characterized  by  a  more  intense  increase  of  iY2 
for  high  values  of  R.  The  geographic  location  of  these  types  of  solar  varia¬ 
tion  of  f  F2  with  R  are  dependent  on  local  time  and  season.  Vasil'yeva 
examined  the  percentage  distribution  of  the  various  types  of  dependencies  of 
f  F2  on  R.  He  subdivided  type  II  into  two  additional  categories.  Type  Ha 
were  cases  where  fQF2  continued  to  grow  at  large  values,  but  more  slowly  than 
at  low  and  medium  solar  activity.  Type  lib  includes  dependencies  where  not 
only  the  growth  of  fQF2  is  limited,  but  also  where  the  critical  frequencies 
decrease  somewhat  with  increasing  solar  activity.  The  percentage  of  cases  of 
each  type  of  dependence  of  fQF2  on  R  were  given  as  follows:  (I)  20.7%,  (II) 
64.8%,  (Ila)  4.2%,  (lib)  3.2%,  and  (III)  7.3%. 

This  improved  version  of  MINIMUF-3.5,  MINIMUF-B,  was  developed  by 
altering  MINIMUF-3.5  so  that  a  version  called  MINIMUF-A  was  developed  for 
which  the  predicted  MUFs  were  not  dependent  on  sunspot  number.  The  equation 

for  f  F2  in  MINIMUF-A  became 

o 

foF2  ■  (Ao  +  A1  'cos  *eff)1/2  <12) 
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MINIMUF-B  was  similar  to  MINIMUF-3.5  except  that  the  expression  for  fQF2 
was  given  by 

foF2  =  f(R)(Ao  +  Aj  ^~K^f)U2  (13) 
where  f(R)  was  a  function  dependent  on  the  sunspot  number  R. 

The  function  f(R)  was  found  by  calculating  the  residuals  between  the 
predicted  MUFs  from  MINIMUF-A  and  the  corresponding  MOFs  from  the  oblique 
sounder  data  base.  The  residuals  were  subdivided  into  subsets  according  to 
sunspot  variation.  The  first  subset  was  for  data  for  which  the  sunspot  number 
used  for  the  calculation  of  the  MUF  varied  from  10  to  20.  In  each  subsequent 
subset  the  sunspot  number  was  incremented  by  10.  The  average  residual  (the 
bias)  is  calculated  in  each  subset  to  help  determine  the  needed  sunspot 
variation.  In  particular,  a  multiplying  constant  A  is  determined  in  each 
sunspot  interval  that  made  the  bias  zero. 

A  polynomial  equation  was  fit  in  a  least  square  sense  to  these  constants. 
For  this  purpose  a  standard  subroutine  for  a  Hewlett-Packard  hand  held 
programmable  calculator  was  used.  The  program  approximated  the  function  f(R) 
by  a  polynomial  of  degree  m,  where  2<m<4.  The  special  Chebyshev  polynomials 
for  discrete  intervals  were  used.  Figure  4  shows  the  results  of  fitting  both 
a  second  order  polynomial  and  a  fourth  order  polynomial.  It  also  shows  the 
multiplier  used  by  MINIMUF-3.5  and  the  multiplier  determined  by  calculating 
the  residuals  between  MUFs  predicted  by  MINIMUF-A  and  MOFs  from  the  observed 
data  base.  First  note  the  oscillatory  nature  of  the  multiplying  constants  A, 
particularly  at  sunspot  numbers  greater  than  100.  The  fourth  order  fit  has  an 
oscillatory  nature  too.  The  second  order  fit  and  MINIMUF-3.5  gives  the  best 
fit  below  sunspot  number  100.  The  second  order  fit  also  has  the  ability  to 
limit  the  growth  of  fQF2  at  sunspot  numbers  above  100. 

A  further  investigation  was  made  to  determine  the  cause  of  the  large  peak 
in  the  observed  multiplying  constants  A.  The  multiplying  constants  are 
determined  as  a  function  of  sunspot  number  for  individual  paths.  It  is  found 
that,  for  two  paths  in  the  oblique  sounder  data  base,  the  fQF2  decreases  with 
increasing  R  for  R  >  (120-140).  This  type  of  dependence  on  sunspot  number 
only  occurred  3.2  percent  of  the  time  (reference  16).  Hence,  the  data  for 


21 


2.25 


2-i 


Figure  4.  Initial  sunspot  number  multiplier  fit. 

these  two  paths,  Winnipeg  to  Resolute  Bay  and  Ottawa  to  The  Hague,  were 
dropped  from  the  oblique  sounder  data  base  used  to  determine  f(R).  Subsequent 
work  on  improving  MINIMUF  in  the  polar  region  showed  these  two  paths  to  be 
particularly  illustrative  of  the  weakness  of  MINIMUF  in  the  polar  regions. 

A  new  second  order  fit  was  made  to  the  multipliers  using  the  adjusted 
oblique  sounder  data  base.  The  multiplier  for  MINIMUF-3.5,  the  new  second 
order  fit,  and  the  data  points  are  shown  in  figure  5.  As  can  be  seen  in 
figure  5,  below  a  sunspot  number  of  about  125  the  second  order  fit  and 
MINIMUF-3.5  differ  only  by  a  small  amount.  But  at  high  sunspot  levels,  as 
occurred  during  July  and  August  1982,  the  second  order  fit  would  provide  a 
limit  on  the  growth  of  the  fQF2  part  of  the  MUF  computation.  Table  8  gives 
the  coefficients  for  the  second  order  fit. 
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SUNSPOT  MULTIPLIER  USED  IN  MININUF 
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Figure  5.  Second  order  sunspot  number  multiplier  fit. 


To  further  improve  the  fit,  the  multiplying  constants  A  were  smoothed  by 
averaging  each  pair  of  data  points  starting  at  the  lowest  sunspot  interval. 
The  resulting  fit  along  with  the  multiplier  from  MINIMUF-3.5  is  presented  in 
figure  6.  Table  8  gives  the  coefficients  for  the  smoothed  second  order  fit. 


The  accuracy  of  both  MINIMUF-3.5  and  the  smoothed  version  of  MINIMUF-B 
was  determined  by  comparing  the  bias  and  the  rms  error  of  the  residuals  from 
both  programs  as  a  function  of  the  sunspot  number.  Overall  the  bias  and  rms 
error  for  MINIMUF-3.5  were  0.52  MHz  low  and  4.33  MHz,  respectively.  For 
MINIMUF-B,  the  bias  and  rms  error  were  0.14  MHz  low  and  4.33  MHz, 
respectively.  To  determine  the  accuracy  of  the  program  as  a  function  of 
sunspot  number,  the  data  were  divided  into  sunspot  number  categories.  Table  9 
gives  the  bias  and  rms  error  for  each  model.  A  negative  number  means  that  the 
model  predicted  high.  As  can  be  seen,  the  most  notable  differences  in  accuracy 
as  a  function  of  the  sunspot  number  are  in  the  bias.  In  the  sunspot  number 
range  10-120,  the  bias  is  always  less  than  1.0  MHz.  However,  at  the  sunspot 
numbers  greater  than  120,  much  higher  biases  are  shown. 
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SUNSPOT  MULTIPLIER  USED  IN  HININUF 


+  SMOOTHED  DATA 


Figure  6.  Second  order  sunspot  number  multiplier  fit  to  smoothed  data. 


,  2 

Table  8.  Coefficients  for  f(R)  =  d0  +  di  R  +  d2R  * 


Version 

d 

0 

dl 

d2 

unsmoothed 

smoothed 

0.966184943 

0.965142201 

0.006147659 

0.006288535 

i 

-0.000014123 

-0.090016288 

Table  9.  Accuracy  of  MINIMUF-3.5  and  the  smoothed  version  of  MINIMUF-B 
as  a  function  sunspot  number. 
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In  the  development  of  MINIMUF-B,  two  paths,  Winnipeg  to  Resolute  Bay  and 
Ottawa  to  The  Hague,  were  deleted  from  the  data  set  to  'obtain  f(R).  With 
these  paths  deleted,  the  overall  accuracy  of  MINIMUF-3.5  was  0.34  MHz  low  with 
an  rms  error  of  4.23.  For  MINIMUF-B,  the  overall  accuracy  was  0.05  MHz  high 
with  an  rms  error  of  4.26  MHz. 

CHOICE  OF  CONTROL  POINTS 

In  MINIMUF-3.5  the  calculated  MUF  is  the  minimum  value  evaluated  at 
specific  "control  points"  along  the  great  circle  propagation  path.  In  the 
original  version  of  MINIMUF-3.5,  the  control  point  locations  for  path  lengths 
greater  than  4000  km  were  located  2000  km  from  either  terminus  (reference  2). 
However,  in  the  next  versions  (references  3  and  18),  the  control  points  for 
path  lengths  greater  than  4000  km  were  located  at  one  quarter,  one  half,  and 
three  quarters  of  the  path  length  along  the  great  circle  path.  In  a  version 
published  in  the  United  Kingdom  (reference  19),  the  control  points  were 
located  one  quarter  and  three  quarters  of  the  path  length  along  the  great 
circle  path.  For  all  versions,  the  control  point  for  path  lengths  less  than 
4000  km  is  located  at  the  path  mid-point.  The  version  of  MINIMUF-3.5  used  to 
determine  the  solar  variation  presented  earlier  in  this  report  was  the 
original  version  with  control  points  located  2000  km  from  either  terminus  for 
path  lengths  greater  than  4000  km. 

This  section  examines  these  three  choices  of  control  points.  '"he 
original  version  of  MINIMUF-3.5  was  modified  so  that  it  calculated  control 
points  at  1/4,  1/2,  and  3/4  of  the  great  circle  path  length;  this  version  was 
called  MINIMUF-D.  Another  version  was  created  that  calculated  control  points 
at  1/4  and  3/4  of  the  great  circle  path  length  only;  this  version  was  called 
MINIMUF-E.  These  two  versions  are  compared  to  the  original  MINIMUF-3.5  as  a 
function  of  path  length  in  tables  10  and  11.  Table  10  gives  the  bias,  and 
table  11  gives  the  rms  error  of  the  models.  From  0-4000  km  only  the  midpoint 
is  used  as  a  control  point  so  there  is  no  difference  between  the  models.  From 
4000  -  6000  km,  the  original  MINIMUF-3.5  had  the  smallest  bias  and  lowest  rms 
error.  However,  from  6000  -  8000  km,  MINIMUF-D  had  the  lowest  bias  of  the 
three  versions  and  a  lower  rms  error  than  MINIMUF-3.5.  The  difference  in  rms 
error  between  MINIMUF-D  and  MINIMUF-E  is  insignificant. 
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Table  10.  Bias  of  MINIMUF  models  with  range  (MHz). 


Range  (km) 

MINIMUF -3. 5 

MINIMUF-D 

MINIMUF -E 

MINIMUF-F 

MINIMUF -I 

0-1000 

0.91 

0.91 

0.91 

0.91 

0.91 

1000-2000 

-2.99 

-2.99 

2.99 

-2.99 

-2.99 

2000-3000 

0.97 

0.97 

0.97 

0.97 

0.75 

3000-4000 

0.20 

0.20 

0.20 

0.20 

0.00 

4000-5000 

0.76 

1.89 

1.87 

0.76 

0.74 

5000-6000 

1.75 

2.34 

2.29 

1.75 

1.58 

6000-7000 

-0.92 

-0.51 

-0.56 

-0.51 

-0.51 

7000-8000 

-0.41 

-0.16 

-0.24 

-0.16 

-0.16 

Table  11.  RMS  error  of  MINIMUF  models  with  range  (MHz). 


Range  (km) 

MINIMUF -3. 5 

MINIMUF-D 

MINIMUF-E 

MINIMUF-F 

MINIMUF -I 

0-1000 

2.34 

2.34 

2.34 

2.34 

2.34 

1000-2000 

4.05 

4.05 

4.05 

4.05 

4.05 

2000-3000 

4.08 

4.08 

4.08 

4.08 

3.99 

3000-4000 

5.14 

5.14 

5.14 

5.14 

5.10 

4000-5000 

3.53 

3.97 

3.96 

3.53 

3.53 

5000-6000 

5.18 

5.41 

5.40 

5.18 

5.07 

6000-7000 

3.83 

3.76 

3.77 

3.76 

3.76 

7000-8000 

4.24 

4.18 

4.17 

4.18 

4.18 

A  fourth  version  of  MINIMUF  was  created  called  MINIMUF-F  that  had  control 
points  at  the  path  midpoint  for  path  lengths  less  than  or  equal  4000  km, 
control  points  located  2000  km  from  either  terminus  for  path  lengths  greater 
than  4000  km,  but  less  than  or  equal  6000  km,  and  control  points  located  at 
1/4,  1/2,  and  3/4  of  the  great  circle  path  length  for  path  lengths  greater 
than  6000  km.  The  results  for  this  version  are  presented  in  tables  10  and  11. 
Obviously  it  has  the  lowest  bias  and  rms  error  of  the  four  versions. 

GEOMAGNETIC  LATITUDE  DEPENDENCE 


To  account  for  critical  frequency  separation  between  ordinary  and 
extraordinary  traces,  it  is  common  in  prediction  methods  to  add  one-half  the 
gyrofrequency  to  the  fQF2  in  determination  of  the  MUF.  However,  in  the  use  of 
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MINIMUF-3.5,  the  determination  of  the  constants  in  the  fitting  process 
includes  implicitly  the  gyrofrequency  for  the  paths  for  Which  the  constants 
were  determined.  For  the  northern  most  path  used  in  the  fitting  process,  the 
gyrofrequency  at  path  midpoint  is  approximately  1.2  MHz.  On  the  avetage,  the 
gyrofrequency  is  approximately  1.0  MHz  in  the  MINIMUF  fitting  process.  At 
high  latitudes  the  gyrofrequency  can  get  as  high  as  1.8.  Hence,  a  bias  can 
exist  in  the  fQF2  as  large  as  0.4  MHz.  For  an  M  factor  of  3.2,  this  causes  a 
bias  in  the  MUF  of  1.28  MHz  low. 


Because  of  this  potential  bias  in  MINIMUF-3.5  at  high  latitudes,  the 
effect  of  adding  one-half  the  gyrofrequency  to  the  fQF2  was  investigated.  Two 
versions  of  MINIMUF  were  produced.  The  first,  MINIMUF-G,  added  one-half  the 
gyrofrequency  to  the  fQF2  at  all  latitudes.  The  second,  MINIMUF-H,  added  one 
half  the  gyrofrequency  to  the  fQF2  at  latitudes  greater  than  55°N  geomagnetic. 
In  each  version  0.5  MHz  was  subtracted  off  to  remove  the  implicit  fitting  of 
the  gyrofrequency  in  MINIMUF-3.5. 


From  Davies  (reference  20),  the  gyrofrequency  for  an  earth  centered 
dipole  field  is  given  by 


fH  =  0.870 


1+3  sin20  )1/Z  (MHz) 


(14) 


where  R  =  earth  radius  (6371  km) 
e 


H_  *  height  of  the  maximum  ionization  of  the  F-layer  at  the  midpoint 
r 


of  the  propagation  path 


0  ®  latitude  of  the  midpoint  of  the  propagation  path  in  magnetic 
coordinates  (radians). 


Substituting  the  approximate  value  for  Rfi  and  300  km  for  R^,  equation  (14) 
becomes 


f„  =  0.75781  1  +  3  sin  0 
H 


)^2  (MHz) 


(15) 
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From  Davies  (reference  20)  the  geomagnetic  latitude  9  is  given  by 

sin  8  *  sin#sin#Q  +  cos#cos#ocos^X  -  \Q ^ 
where  #  =  latitude  of  the  midpoint  of  the  propagation  path  (radians) 

X  =  longitude  of  the  midpoint  of  the  propagation  path  (radians) 

A  *  latitude  of  the  North  magnetic  pole 
o 

(1.3666  radians,  North  or  78.3°N) 

X  =  longitude  of  the  North  magnetic  pole 
(1.2043  radians,  West  or  69°W). 

Substituting  the  values  of  4»o  and  \q  into  equation  (17)  provides 

sine  =  0.9792  sin*  +  0.2028  cos#  cosU  -  1.2043)  .  (17) 

The  accuracy  of  MINIMUF-G,  MINIMUF-H,  and  MINIMUF-3.5  are  compared  to 
observed  MOFs  on  the  39  paths.  The  latitude  of  the  control  points  are  divided 
into  five  latitude  regions:  transequator ial  (Th),  low  latitude  (LO),  mid¬ 
latitude  (M),  high  latitude  (H),  and  transauroral  (TA).  Table  12  shows  the 
path  categories  for  the  additional  paths  added  to  the  data  base  since  the  1981 
comparison  (reference  4).  The  path  characteristics  for  the  first  25  paths  are 
given  in  Table  10  of  reference  4.  Tables  13  and  14  show  the  bias  and  rms 
error  of  these  models,  respectively.  Only  at  high  and  transauroral  latitudes 
is  there  an  impro’  nt  in  predicted  MUF  by  using  one-half  of  the  gyro- 
frequency  added  to  the  fQF2. 

MINIMUF-I  included  in  tables  10,  11,  13,  and  14  is  a  combination  MINIMUF- 
F,  optimum  control  point  version,  with  MINIMUF-H,  half  the  gyrofrequency  added 
to  f  F2  at  high  latitudes.  Note  that  MINIMUF-I  is  less  accurate  in  trans¬ 
equator  ial  and  low  latitude  regions  even  though  MINIMUF-H  is  better  there. 
This  occurs  because  on  three  very  long  paths  the  additional  control  point  at 
path  midpoint  causes  an  additional  bias  in  the  result.  The  accuracy  of  the 
prediction  at  the  control  points  at  path  midpoint  is  probably  less  accurate 
than  those  control  points  at  higher  latitudes  nearer  the  path  terminals. 
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Table  12.  Additional  path  characteristics. 


No.  Transmission  Path  Orientation 

26  Puerto  Rico  to  Maynard,  MS  N-S 

27  Thule,  Greenland  to  Stockbridge,  NY  N-S 

28  Andoya,  Norway  to  Maynard,  MS  E-W 


29  Bangkok,  Thailand  to  Chantaburi,  Thailand  Other 

30  Ottawa,  Canada  to  The  Hague,  Netherlands  Other 

31  Winnipeg,  Canada  to  Resolute  Bay,  Canada  N-S 

32  Ottawa,  Canada  to  Resolute  Bay,  Canada  N-S 

33  Okinawa  to  St.  Kilda,  Australia  N-S 

34  Okinawa  to  Townsville,  Australia  Other 

35  Yamagawa,  Japan  to  St.  Kilda,  Australia  N-S 

36  Yamagawa,  Japan  to  Townsville,  Australia  Other 

37  Monrovia,  Liberia  to  Rota,  Spain  N-S 

38  Monrovia,  Liberia  to  Fort  Lamy,  Chad  Other 

39  Tripoli,  Libya  to  Accra,  Ghana  Other 

TE  ■  Transequatorial  E-W  * 

LO  “  Low  latitude 

M  =  Mid-latitude 
H  *  High  latitude 
TA  =  Transauroral 


Latitude  of 
Control  Points 


M 
TA 
TA 
LO 
H 
TA 
TA 
TE 
TE 
TE 
TE 
LO 
LO 
LO 

East/West 


Region 


B 

C 

C 

A 

C 

A 

A 

C 

C 

C 

B 

A 

A 

A 


N-S  =  North/South 
A  =  Continental 
B  ■  Ocean 

C  =  Combined  land /ocean 


Table  13.  Bias  of  MINIMUF  models  with  geomagnetic  latitude  region  (MHz). 


Latitude  Region 

MINIMUF- 3. 5 

MINIMUF-G 

MINIMUF-H 

MINIMUF- I 

TE 

0.31 

0.63 

0.31 

0.43 

LO 

0.87 

1.13 

0.87 

1.12 

M 

-0.24 

-0.38 

-0.24 

-0.19 

H 

1.98 

1.54 

1.60 

1.60 

TA 

1.73 

1.14 

1.14 

1.14 
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Table  14.  RMS  error  of  MINIMUF  models  with  geomagnetic  latitude  region  (MHz). 


Latitude  Region 

MINIMUF -3. 5 

MINIMUF-G 

MINIMUF -H 

MINIMUF -I 

TE 

4.85 

4.91 

4.85 

4.80 

LO 

4.79 

4.83 

4.79 

4.83 

M 

3.47 

3.48 

3.47 

3.44 

H 

4.12 

3.93 

3.96 

3.96 

TA 

5.41 

5.25 

5.25 

5.25 

Table  15  shows  the  progression  of  improvement  in  the  models  with  each 
succeeding  change.  MINIMUF-I  is  the  version  of  MINIMUF  for  which  further 
development  was  done. 


Table  15.  Overall  comparison  of  MINIMUF  models. 


Program 

Biaz 

MHz 

rms  error 

MHz 

MINIMUF -D 

0.98 

4.46 

MINIMUF -E 

0.95 

4.46 

MINIMUF -F 

0.60 

4.32 

MINIMUF-G 

0.48 

4.32 

MINIMUF -H 

0.42 

4.29 

MINIMUF -I 

0.51 

4.28 

MINIMUF -3. 5 

0.51 

4.33 

CRITICAL  FREQUENCY  MODEL 


In  this  section  the  development  of  an  improved  fQF2  portion  of  MINIMUF  is 

presented.  This  involves  the  fitting  of  the  expression  for  fQF2  portion  to 

vertical  incidence  f  F2  data. 

o 

Equation  (3)  presented  earlier  shows  that  the  only  sunspot  number 

dependence  in  the  MINIMUF  3.5  model  is  through  the  linear  multiplying  factor 

1  +  R/R  .  The  value  of  R  =  250  is  chosen  based  on  the  limited  data  set  used 
o  o 

to  determine  the  contents  in  the  model  as  described  above.  This  data  set  had 
a  maximum  sunspot  number  of  85.  Hence,  it  is  perhaps  not  surprising  that 
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the  model  proved  less  than  adequate  for  the  range  of  sunspot  numbers 
encountered  during  July/ August  1982  test. 

The  original  MINIMUF  3.5  algorithm  is  derived  using  the  12  month  running 
mean  sunspot  number.  In  practice  however,  all  types  of  sunspot  numbers  have 
been  used  as  input  to  MINIMUF- 3. 5  for  prediction  purposes.  During  the  July/ 
August  1982  test,  for  example,  daily  observed  values  of  sunspot  numbers  were 
input  to  MINIMUF  with  much  less  than  satisfactory  results.  Since  this 
procedure  is  followed  in  many  instances,  especially  when  MINIMUF  is  used  as  a 
real-time  frequency  selection  aid,  we  have  changed  the  sunsp.  t  representation 
from  the  smoothed  12  month  running  mean  sunspot  number  to  the  monthly  median 
values.  This  more  closely  agrees  with  what  is  used  in  practice.  However,  one 
can  still  expect  significant  errors  when  daily  sunspot  number  values  are 
input.  This  is  especially  true  during  times  of  ver'>  active  and  disturbed 
solar  conditions. 

In  the  development  that  follows,  the  basic  functional  form  of  equation 

(3)  is  retained  with  "he  attempt  to  develop  a  new  model  altogether.  The 

predictive  capability  of  the  current  model  is  quite  good  for  low  to  mid-range 

values  of  sunspot  number,  and  our  intent  is  to  improve  the  accuracy  for  the 

high  and  very  high  sunspot  range  without  compromising  its  success  in  these 

other  ranges.  With  this  criterion  established,  we  began  an  investigation  of 

the  values  of  A  and  A.  as  they  existed  in  MINIMUF-3.5. 
o  1 

We  began  by  graphically  comparing  the  diurnal  predictions  of  equation  (3) 
to  vertical  incidence  fQF2  data  gathered  at  several  representative  sites 
described  earlier  for  various  sunspot  number  values.  These  results  showed 
that,  in  general,  the  bias  value,  Aq>  in  equation  (3)  is  not  a  strong  function 
of  sunspot  number.  The  value  of  A^,  however,  allowed  too  much  diurnal 
amplitude  variation  at  very  high  sunspot  number  values.  Figure  7  illustrates 
the  effect  of  the  value  of  A^  being  too  large.  In  the  case  presented  (Moscow, 
August  1970,  SSN  =  194),  the  maximum  predicted  value  is  about  2.5  MHz  high. 
An  obvious  approach  then  is  to  remove  the  linear  function  of  sunspot  number  in 
equation  (3),  and  let  A^  become  the  only  function  of  sunspot  number  in  the  new 
fQF2  expressions.  This  function  of  sunspot  number  replaced  that  derived  for 
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MINIMUF -B  presented  earlier.  The  fQF2  data  set  was  then  partitioned  into  24 
segments,  eac*'  containing  a  range  of  10  in  sunspot  number,  from  the  lowest 
value  in  the  data  base,  approximately  1,  to  the  highest  of  approximately  245. 
For  each  of  these  segments,  values  of  [(fQF2^)  -  Aq]  are  compared  to  values 
of  /cos  Xgff»  where  cos  xeff  is  the  value  returned  from  MINIMUF  3.5  for  the 
particular  time  and  location  corresponding  to  each  data  point  f oF2^ . 

is 
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Figure  7.  Example  of  the  comparison  of  MINIMUF-3.5  fQF2  with  vertical 
iorosonde  data. 


COMPARISON  OF  MINIMUF  3.S  F0F2  WITH  SOUNDER  DATA 
MOSCOW i  AUGUST  1970  SSN-194 


The  result  of  the  comparison,  using  the  DASCR3  statistical  data  compari¬ 
son  program,  is  a  set  of  24  mean  square  estimated  values,  one  for  sunspot 
number  segment,  of  the  multiplier  A^,  where 


(f  F2.)2 
o  a 


/cos  x 


eff ' 


(18) 


A  linear  least  squares  fit  is  then  derived  for  this  set  of  24  data 
points,  resulting  in  the  equation 


A^SSN)  =  0.814R  +  22.23  .  (19) 

Because  A^  appears  under  the  square  root  sign  in  equation  (3),  the  resulting 
expression  for  fQ12  is  nonlinear  with  sunspot  number.  Figure  8  shows  the 
results  using  the  new  expression  versus  that  obtained  from  MINIMUF- 3. 5( I )  for 
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a  cos  xeff  ■  0.5.  Note  that  the  new  expression  provides  a  saturation  effect, 
albeit  a  slow  one,  in  the  behavior  of  the  critical  frequency  as  a  function  of 
sunspot  number. 


F 

o 

F 

2 

I 

N 

H 

H 

Z 


COMPARISON  OF  MINIMUF  3.3  F0F2  WITH  HEM  MODEL  VERSUS  SUNSPOT  HUMBER 
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Figure  8.  Comparison  of  MINIMUF-3.5  fQF2  with  new  model  versus  sunspot 
number  f-*r  cos  xeff  =  0.5. 


A  comparison  of  the  predictions  of  this  improved  fQF2  expression  versus 
data,  for  each  of  the  four  ranges  of  sunspot  number  and  for  the  entire  data 
base  is  given  in  table  16.  In  this  table  the  first  value  shows  the  results  of 
a  comparison  of  the  MINIMUF-3.5(I)  fQF2  expression,  equation  (3),  in  each  of 
these  sunspot  number  ranges.  The  average  difference  (the  bias),  the  rms 
error,  and  the  correlation  coefficient  between  the  models  and  the  observed 
data  are  given. 

The  results  of  the  comparison  show  that  the  changes  to  the  critical 
frequency  portion  of  the  MINIMUF-3.5  algorithm  significantly  improves  the 
overall  accuracy  of  the  prediction.  Note  particularly  the  change  in  the  cor¬ 
relation  coefficient  between  the  models  and  the  f  F2  values.  This  shows 
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clearly  that  the  new  expression  is  a  much  better  representation  of  fQF2.  The 
correlation,  however,  does  drop  somewhat  with  increasing  sunspot  number. 

Table  16.  Statistical  Comparison  of  Accuracy  of  fQF2 
Predictions  of  MINIMUF-3.5  and  MINIMUF-85. 


MINIMUF-3.5  /  MINIMUF-85 


SSN  RANGE 

AVERAGE  DIFF 
(0BS-PRED) 

RMS  VALUE,  DIFF 
(MHZ) 

C0RR.  COEFF. 

Low  ( <30 ) 

-0.97/0.44 

2.51/1.56 

0.24/0.65 

Med  (31-100) 

-0.71/0.15 

2.92/1.86 

0.21/0.63 

High  (101-150) 

-0.43/0.33 

3.66/2.55 

0.14/0.58 

Very  high  (>150) 

-0.94/0.15 

4.29/3.25 

0.05/0.50 

Entire  data  base 

-0.77/0.26 

3.37/2.35 

0.39/0.69 

These  improvements  in  the  f  F2  model  will  be  reflected  most  dramatically 
in  the  short  path  length  MUF  prediction  of  the  new  model  where  the  fQF2 
dominates  over  the  M-f actor.  We  will  also  see  improvements  in  ray  trace 
applications  where  accurate  critical  frequencies  are  important  in  predicting 
mode  structure  when  operating  close  to  the  MUF. 

M- FACTOR  MODEL 

In  the  previous  section,  a  new  f  qF2  representation  for  MINIMUF  is 
determined  by  fitting  equation  (3)  against  fQF2  vertical  incidence  sounder 
data.  Having  done  this,  it  was  necessary  to  make  adjustments  to  the  M-f actor 
model  so  that  the  resulting  predicted  MUF  is  in  fact  a  MUF.  This  is  done  by 
comparing  the  resulting  predicted  MUF  to  the  observed  MOF  data  base.  In  so 
doing,  sunspot  number,  seasonal,  and  diurnal  dependencies  are  incorporated 
into  the  M- factor  model. 

SUNSPOT  NUMBER  DEPENDENCE 

As  shown  in  figure  8,  the  changes  in  the  fQF2  algorithm  described  above 
provide  a  "saturation"  effect,  albeit  a  slow  one,  in  the  behavior  of  the 
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critical  frequency  as  a  function  of  the  sunspot  number.  This  is  a  frequently 
observed  effect  in  the  dynamics  of  the  critical  frequency*  (references  15,  16, 
17).  It  is  also  known  that  the  obliquity  factor,  hereafter  called  the 
M-factor,  shows  an  inverse  dependence  on  sunspot  number  (references  21  and 
22).  The  combination  of  these  two  behaviors,  combined  in  the  product  of 
equation  1,  produces  a  natural  and  relatively  fast  cut-off  in  the  use  of  the 
MUF  as  a  function  of  the  sunspot  number.  It  is  felt  that  the  poor  performance 
of  MINIMUF-3.5  during  periods  of  very  high  sunspot  number  can,  at  least 
partially,  be  attributed  to  the  lack  of  sunspot  number  dependence  in  the 
M-factor. 


In  order  to  determine  the  sunspot  number  dependence  of  the  M-factor,  we 
followed  a  similar  procedure  to  that  described  above  for  the  critical 
frequency.  That  is,  the  MOF  data  base  described  earlier,  which  consisted  of 
hourly  median  values  of  MOF  scaled  from  oblique  ionograms  part  of  which 
represented  a  maximum  monthly  median  sunspot  number  of  approximately  145,  was 
subdivided  into  15  segments,  each  defined  by  a  range  of  10  in  sunspot  number. 
Using  the  DASCR3  statistical  program,  we  compared  prediction  using  equation  1, 
using  the  improved  fQF2  described  above  instead  of  the  original  fQF2,  with 
actual  data  values.  This  resulted  in  a  set  of  15  multipliers,  A^,  defined  by 


MOF,  =  A„(SSN)  x  M  x  /A  +  A.(SSN)  *  /cos  x 
a  L  o  1  err 

Again  a  linear  least  square  fit  is  made  to  the  A 2  data  points. 

resulted  in  the  function  for  A^, 

A2(SSN)  =  1.3022  -  0.00156  R 


(20) 

This 


(21) 


which,  as  expected,  shows  the  raonotonically  decreasing  behavior  as  a  function 
of  the  sunspot  number.  Figure  9  shows  a  comparison  of  the  derived  multiplying 
factors  and  the  linear  fit  to  the  factors.  Figure  10  shows  a  comparison  of 
this  new  MUF  expression  versus  the  original  as  a  function  of  the  sunspot 
number  for  a  3000  km  path  with  cos  xe^£  =0.5. 


It  should  be  noted  that  figure  10  shows  saturation  at  a  MUF  value  of 
approximately  11  MHz  for  this  case.  This  relatively  low  value  is  due  to  a 
lack  of  a  very  high  sunspot  number  data  in  the  MOF  data  base.  It  would  be 
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COMPARISON  OF  DERIVED  MULTIPYING  FACTORS  AND 
LINEAR  FIT  TO  FACTORS 


Figure  9.  Comparison  of  derived  sunspot  multiplying  factors  and  a  linear 
fit  to  the  factors. 


COMPARISON  OF  HINIHUF  3.3  WITH  HEW  MODEL  VERSUS  SUNSPOT  NUMBER 
3000  KM  PATH  WITH  COS  XL- >.3 


Figure  10.  Comparison  of  MINIMUF-3.5  with  new  model  versus  sunspot 
number,  3000  km  path  with  cos  xeff  =  0.5. 
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desirable  to  obtain  more  high  sunspot  number  data  (>  150).  The  inclusion  of 
this  data  and  a  repeat  of  the  process  described  above,  should  lead  to  a  higher 
cut-off  value  for  the  MUF. 

The  important  point  to  be  stressed  is  that  the  technique  described  does, 
in  fact,  lead  to  a  cut-off  in  the  MUF  with  increasing  the  sunspot  number. 
Th- f.  is  in  line  with  observation  which  shows  that,  unlike  the  linear  behavior 
described  by  equation  (1),  the  MUF  does  saturate  at  high  sunspot  number 
values. 

In  table  17  we  show  the  effect  of  these  changes,  both  in  the  fQF2  and  the 
M-f actor,  on  the  prediction  of  the  MUF  for  the  limited  range  of  sunspot 
numbers  in  the  MOF  data  base.  Note  that  these  changes  have  not  significantly 
altered  the  overall  predictive  capability  of  the  algorithm  for  these  smaller 
sunspot  number  values,  as  intended.  The  increase  in  rms  error  in  the  new 
model  might  be  due  to  the  use  of  a  linear  fit  to  the  multiplying  factors 
rather  than  a  second  order  fit. 


Table  17.  Comparison  of  accuracy  of  MINIMUF-3.5  with  intermediate 
version  of  new  model  with  SSN  dependence  in  fQF2  and  M- 
f actor . 


Avg.  Diff. 

RMS  Diff. 

Corr.  Coef. 

(MHz) 

(MHz) 

MINIMUF-3.5 

0.52 

4.28 

0.85 

INTERMEDIATE 

-0.11 

4.52 

0.85 

SEASONAL  DEPENDENCE 

It  is  well  known  that  the  MOF  is  a  seasonally  dependent  parameter. 
During  the  winter  months,  for  example,  increased  electron  density  levels  in 
the  ionosphere  generally  cause  a  lowering  of  the  height  of  maximum  electron 
density.  This  "winter  anomally"  allows  higher  frequencies  to  propagate  on  a 
given  transmission  path.  Consequently,  we  normally  see  larger  MOFs  in  the 
winter  than  in  the  summer  or  equinox  months.  To  include  these  effects  in  the 
model,  we  followed  a  similar  procedure  to  that  described  above  to  derive  a 
seasonal  (monthly)  dependence  in  the  M-factor  part  of  the  algorithm. 
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The  MOF  data  base  is  again  segmented,  this  time  by  month,  into  12 
separate  data  bases.  For  each  of  these  data  bases,  we  compare  prediction  from 
equation  (1),  with  the  above  changes  incorporated,  to  the  measured  data.  That 
is,  we  compared  the  predictions  of  the  expression 


A-(month)  *  A-(SSN)  *  M  *  f  F2 
j  L  o 

to  corresponding  values  of  MOF^  from  the  data  base  using  DASCR3.  This 

resulted  in  12  values  of  the  multiplier  A_,  one  for  each  month.  These  values 

th  ^ 

were  then  fit  with  a  6  order  Fourier  series  given  by 


A^(month)  =  0.9925  +  O.OllsinM  +  0.087cosM 

-  0.043sin2M  +  0.003cos2M 

-  0.013sin3M  -  0.022cos3M 
+  0.003sin4M  +  0.005sin5M 
+  0.018cos6M 


where 


M 


2tt  month 
12 


(22) 


In  figure  11  we  show  the  derived  multiplier  and  the  Fourier  fit  to  the 
multipliers.  Note  that  there  is  a  small  increase,  compared  to  the  summer 
months,  evident  in  these  multiplers,  thus  giving  a  small  increase  to  the  MUF 
values  during  the  winter  period.  In  table  18  we  show  a  comparison  of  the 
predictions  of  this  intermediate  algorithm,  with  seasonal  and  the  sunspot 
number  dependence  in  the  M- factor  and  improved  fQF2,  against  those  of 
MINIMUF-3.5  for  the  entire  MOF  data  base.  Again  we  see  that,  for  the  ranges 
of  sunspot  number  and  seasons  contained  in  the  MOF  data  base,  the  results  of 
the  above  changes  have  not  significantly  altered  the  predictive  capability  of 
the  original  algorithm.  We  feel  that  with  the  seasonal  dependence  in  the  M- 
f actor,  we  have  a  more  versatile  and  accurate  prediction  program  which 
reflects  more  of  the  dependencies  we  see  in  the  MOF  parameter.  What's  more, 
we  also  have  a  more  accurate  M-factor  which  can  be  used  for  other  purposes 
within  the  PROPHET  program. 
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MONTHLY  MULTIPLYING  FACTOR  AND  FOURIER  FIT 
M  FACTOR 
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Figure  1 1 .  Monthly  multiplying  factor  and  Fourier  fit  to  M-factor. 


Table  18.  Comparison  of  accuracy  of  MINIMUF-3.5  with  intermediate 
version  with  seasonal  dependence  in  M  factor. 


Avg.  Diff. 

RMS  Diff. 

Corr.  Coef. 

(MHz) 

(MHz) 

MINIMUF-3.5 

0.52 

4.28 

0.85 

INTERMEDIATE 

0.02 

4.15 

0.85 

TIME  DEPENDENCE 

Finally,  in  order  to  have  a  diurnally  varying  M-factor  which  would 
reflect  diurnal  changes  in  the  height  of  the  F2  layer,  we  used  the  procedure 
described  above  to  determine  the  time  dependence  of  the  M-factor. 

In  order  to  separate  night  from  day  in  the  paths  that  make  up  the  MOF 
data  base,  it  is  of  course  necessary  to  use  local  time  at  the  control  points 
for  this  p.  sdure.  Then,  following  the  method  of  the  other  cases  described 
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above,  the  data  base  is  segmented  into  two  parts,  a  seveu  hour  "day"  and  a 
"night"  part.  A  comparison  of  the  predictions  of  the  expression 

A,  (time)  *  A,  *  A»  x  M  *  f  F2 
4  3  2  o 

with  the  day  part  of  the  data  base  led  to  a  set  of  multipliers,  A^,  which 
could  be  adequately  fit  with  a  linear  function  of  local  time. 

(DAY)  A^(time)  =  1.11  -  0.01  r  (23) 

For  the  night  part  of  the  procedure,  some  problems  arose.  The  complica¬ 
tion  is  that  length  of  day  (and  night)  differs  greatly  among  the  paths  in  the 

data  base.  This  leads  to  a  large  amount  of  fluctuation  in  the  day/night 

transition  time  multipliers.  It  proved  necessary  to  introduce  a  new  time 
coordinate,  hours  after  sunset,  in  order  to  adequately  fit  the  nighttime 

multipliers.  This  fit  is  accomplished  with  a  6t^1  order  Fourier  series, 

(NIGHT)  A4(time)  =  1.0195 

-0.06sin2t  -  0.037cos2t 
+0.018sin4t  -  0.003cos4t 
+0.025sin65  +  0.018cos6t 
+0.007sin  8t  -  0.005cos8t 
+0.006sinl0t  +  0.017cosl0t 
-0.009sinl2t  -  0.004cosl2t 

where  t  =  t,  ,  -  t 

local  sunset 

Table  19  shows  a  final  comparison  of  the  new  version  of  the  MINIMUF 
algorithm,  MINIMUF-85,  to  the  prediction  of  MINIMUF-3.5  against  the  MOF  data 
base.  We  note  that  there  is  a  significant  improvement  in  the  overall  bias  of 
the  model  and  a  relatively  small  improvement  in  the  RMS  value  and  correlation 
coefficient.  Overall,  then,  the  predictive  capability  of  the  model  has 
improved  even  over  the  relatively  restricted  set  of  data  in  the  MOF  data  base. 

Table  3  9.  Statistical  comparison  of  MINIMUF-3.5  and  MINIMUF-85 
predictions  of  MUF  against  MOF  data  base. 


Avg.  Diff. 

RMS  Diff. 

Corr.  Coef. 

(MHz) 

(MHz) 

MINIMUF-3.5 

0.52 

4.28 

0.85 

MINIMUF-85 

0. 16 

4.19 

0.86 
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POLAR  REGION  CRITICAL  FREQUENCY  MODEL 


The  physical  basis  for  MINIMUF-3.5(I)  is  that  all  variations  of  the  free 
electron  density  in  the  ionosphere  are  driven  by  solar  zenith  angle.  MINIMUF 
is  developed  by  starting  with  this  assumption,  developing  a  corresponding 
mathematical  model,  and  then  fitting  associated  free  parameters  to  a  set  of 
oblique  sounder  data.  It  is  not  then  surprising  that  MINIMUF-3.5(I)  encounters 
problems  when  predicting  MUFs  at  polar  latitudes.  This  weakness  is  rooted  in 
two  factors.  First,  an  important  contribution  to  ionization  at  high  latitudes 
is  made  by  particle  precipitation.  Solar  photon  induced  dissociation  is  not 
the  only  significant  source  of  ionospheric  free  electrons.  Second,  data  from 
radio  circuits  with  control  points  in  the  polar  regions  are  not  included  in 
the  parameter  fitting  process  in  the  development  of  MINIMUF.  The  control 
point  with  the  highest  geomagnetic  latitude  included  in  the  determination  of 
MINIMUF  parameters  is  the  midpoint  between  Toulouse,  France,  and  Keflavik, 
Iceland.  The  value  of  the  geomagnetic  latitude  at  this  control  point  is 
58.1°N,  which  is  too  low  to  include  the  effects  of  particle  precipitation. 

The  fact  that  high  latitude  ionospheric  behavior  differs  sharply  from 
that  at  lower  latitudes  has  been  recognized  for  some  time.  This  sharp 
difference  requires  the  introduction  of  new  routines  into  MINIMUF  specifically 
tailored  to  model  the  behavior  of  the  MOFs  at  high  latitudes.  By  merging  the 
specific  polar  model  with  the  model  of  lower  latitude  behavior  already  present 
in  MINIMUF,  significant  improvements  in  the  prediction  of  high  latitude  MOFs 
can  be  realized. 

THE  POLAR  MODEL 

The  Chiu  Polar  Model  (reference  23)  for  the  F2  layer  is  developed  as  one 
component  of  a  global  phenomenological  model  of  ionospheric  electron  density. 
The  basis  of  the  model  is  an  analysis  by  Arima  and  Gonezawa  (reference  24)  of 
variations  of  electron  density  into  seasonal,  nonseasonal  annual,  and  non- 
seasonal  semiannual  categories.  The  first  version  of  a  global  model  as 
developed  by  Ching  and  Chiu  (reference  25)  separates  the  global  variations 
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into  polar  and  nonpolar  regimes.  The  polar  and  nonpolar  functions  describing 
each  regime  are  welded  together  by  means  of  a  folding  function. 

The  folding  function  determines  when  polar  effects  (particle  precipita¬ 
tion)  become  dominant.  It  is  a  function  of  geomagnetic  latitude  and  the 
sunspot  number.  The  folding  function  makes  a  fairly  abrupt  transition  from  0 
to  1  between  geomagnetic  latitudes  of  60°  to  75°.  Figure  12  is  a  plot  of  the 
folding  function  for  a  sunspot  number  of  zero.  When  the  folding  function  is 
near  one,  particle  precipitation  effects  are  supposed  to  dominate.  When  the 
folding  function  is  near  zero,  solar  zenith  angle  is  the  major  factor  in 
causing  ionization.  In  between  there  exists  a  fairly  narrow  transition  region 
where  both  sources  of  free  electrons  are  significant. 


Figure  1 2.  The  folding  function  for  monthly  smoothed  sunspot  number  -  0. 

The  Chiu  Polar  Model  represents  a  refinement  of  the  Ching-Chiu  model. 
Both  models  include  diurnal,  seasonal,  and  nonseasonal  annual  dependence.  The 
major  improvement  in  the  Chiu  model  involves  generating  distinct  models  for 
the  North  and  South  polar  regions.  This  change  recognizes  some  of  the 
peculiar  longitudinal  dependencies  which  characterize  the  South  polar  region. 
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The  Chiu  Polar  Model  is  based  on  vertical  sounder  data  gathered  from  18 
different  polar  stations  over  an  entire  solar  cycle.  The  model  predicts  large 
scale  variations  in  ionospheric  electron  density  and  thus  may  be  directly 
applied  to  the  task  of  predicting  fQF2. 

PROCEDURE  FOR  INCORPORATING  THE 
CHIU  POLAR  MODEL  INTO  M3NIMUF 


The  Chiu  Polar  Model  predicts  the  value  of  electron  density.  It  is  the 

electron  densities  from  the  two  regimes  that  the  folding  function  is  designed 

to  properly  combine.  In  folding  a  polar  function  into  MINIMUF,  it  is  necessary 

to  isolate  that  portion  of  MINIMUF  which  calculates  fQF2  and  then  converts  the 

value  of  f  F2  into  electron  density.  MINIMUF' s  f  F2  is  found  by  dividing  its 
o  o 

calculated  value  for  MUF  at  each  control  point  by  the  range -dependent  portion 
of  the  M  factor.  Note  that  the  G  factors  which  are  empiricaly  latitude  depen¬ 
dent  adjustments  to  MINIMUF  remain  in  the  value  to  fQF2  produced  by  MINIMUF. 

f  F2  is  then  converted  to  electron  density  by  use  of  the  equation  (25). 
o 

f  F2  (MHz)  *  2.85  N1^2(electrons/cm3)  .  (25) 

o 

The  electron  density  from  MINIMUF  is  multiplied  by  a  factor  of  one  minus 
the  folding  function  and  then  added  to  the  product  of  the  folding  function  and 
the  Chiu  Polar  Model  electron  density  as  shown  in  equation  (26). 


Ntotal  *  *1_f  ^MINIMUF  +  f  Npolar 


(26) 


The  total  electron  density  at  the  control  point  is  then  converted  back  to  an 
fQF2  at  the  control  point  by  using  equation  (25).  Finally  the  MUF  is  obtained 
by  multiplying  the  value  of  fQF2  by  the  range -dependent  portion  of  the 
M-factor. 


PROCEDURES  USED  FOR  TESTING 
MODEL  ACCURACY  AT  POLAR  LATITUDES 

One  of  the  difficulties  associated  with  building  a  polar  model  is  the 
lack  of  adequate  oblique  sounder  data  with  control  points  located  where  polar 


43 


effects  are  important.  The  only  type  of  data  available  from  south  polar 
regions  is  in  the  form  of  fQF2  soundings,  the  study  of  which  formed  the  basis 
of  the  Chiu  Polar  Model  to  begin  with.  Consequently,  no  study  of  the  accuracy 
of  the  south  polar  model  was  attempted. 

For  checking  the  effect  of  introducing  the  north  polar  model  into  the 
calculation  of  the  MUF,  five  different  paths  comprising  a  total  of  52  path 
months  of  MOFs  formed  the  data  base.  In  these  paths  the  value  of  the  folding 
function  at  the  control  points  varied  from  0.20  to  1.0.  For  paths  1  and  2, 
monthly  average  MOFs  are  taken  at  four  hour  intervals  when  available.  For 
paths  3,  4,  and  5  the  sampling  rate  of  the  monthly  average  rate  is  set  at  two 
hour  intervals.  The  data  on  the  specific  paths  sampled  and  the  number  of 
sample  points  are  presented  in  table  20.  The  first  two  paths  are  additional 
paths  not  contained  in  the  39  path  oblique  sounder  MOF  data  base.  The  MOFs 
are  compared  with  the  MUFs  predicted  by  the  following  three  algorithms:  (1) 
MINIMUF-3.5(I);  (2)  MINIMUF-3 . 5(1)  combined  with  the  Chiu  Polar  Model  via  the 
folding  function;  and  (3)  MINIMUF-85  combined  with  the  Chiu  Polar  Model  via 
the  folding  function. 

The  errors  are  computed  using  the  following  formulas: 

N 

Bias  =  £  [  (MUFpvedicted  i  '  M0Fobserved  i) 
i=l 

(  N  2 

RMS  error  =-  H  £  (MUFpredicted  .  -  M0Fobserved  ,)2 

'  i=l 

where  N  =  the  number  of  sample  points.  A  positive  bias  for  the  polar  model 
error  analysis  means  the  model  predicts  high. 

RESULTS  OF  COMPARISON  OF  POLAR  MODEL 

A  review  of  table  20  makes  clear  the  fact  that  MINIMUF-3. 5(1)  contains  an 
inadequate  model  of  the  polar  ionosphere.  This  inadequacy  may  be  traceable  to 
the  fact  that  at  polar  latitudes,  particle  precipitation  becomes  the  dominant 


Table  20. 

Prediction  errors  in  various  versions  of  MINIMUF  over  polar  paths 
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source  of  ionization  while  MINIMUF-3.5(l)  hypothesizes  that  ionization  is 
driven  solely  by  solar  zenith  angle.  It  would  then  be  expected  that 
MINIMUF-3.5(I)  would  produce  its  largest  errors  during  periods  when  particle 
precipitation  effects  are  most  significant. 

This  is  the  case  illustrated  in  the  data  displays  of  figures  13  through 
18.  During  polar  winter,  minimal  direct  solar  illumination  is  available. 
MINIMUF-3.5(I)  therefore  predicts  relatively  low  values  for  MUFs  influenced  by 
high  latitude  control  points.  However,  particle  precipitation  effects 
continue  unabated  during  polar  winter.  Thus  MINIMUF-3.5(I)  should  consistently 
underestimate  the  wintertime  MUFs  along  polar  paths  to  a  greater  degree  than 
during  other  seasons.  As  shown  in  figure  13,  the  winter  errors  are  almost 
double  the  error  produced  during  the  summer.  Particle  precipitation  effect 
also  increase  with  the  sunspot  activity.  Therefore,  the  error  in  MINIMUF 
3.5(1)  should  be  positively  correlated  with  the  sunspot  number.  This 
hypothesis  is  confirmed  by  the  results  displayed  in  figure  16.  The  error  in 
MINIMUF-3.5(I)  rises  dramatically  as  sunspot  number  increases  above  50. 

Figures  14,  15,  17,  and  18  illustrate  that  both  MINIMUF-3.5(I)  enhanced 
by  the  Chiu  Polar  Model  and  MINIMUF-85  enhanced  by  the  Chiu  Polar  Model  also 
produce  their  largest  errors  in  winter  and  at  a  high  sunspot  number.  Part  of 
this  error  may  be  traced  to  the  fact  that  folding  function  for  all  sampled 
paths  includes  some  fraction  of  MINIMUF-3.5(I)  which  results  in  folding  in  a 
fractional  amount  of  the  error  from  MINIMUF-3.5(I) .  Figure  18  indicates  that 
the  improved  M  factor  present  in  MINIMUF-85  improves  the  performance  of  MUF 
prediction  with  variances  in  the  sunspot  number. 

MINIMUF-3.5(l)'s  estimates  of  MUF  are  an  overall  average  of  about  7  MHz 
low  over  the  polar  paths  sampled.  This  underestimate  is  improved  to  about  2 
MHz  by  adding  the  Chiu  Polar  Model.  By  enhancing  MINIMUF-85  with  the  Polat 
Model,  the  overall  average  bias  is  reduced  to  about  0.5  MHz.  These  results 
again  support  the  contention  that  the  particle  precipitation  effects  ignored 
in  MINIMUF-3. 5(1)  are  significant. 
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Figure  1 5.  MINIMUF-85  folded  with  polar  model  rms  error  as  a  function  of 
month  for  five  selected  polar  paths. 
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Figure  16.  MINIMUF-3.5(I)  rms  error  as  a  function  of  sunspot  number  for 
five  selected  polar  paths. 
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Figure  17.  MINIMUF-3.5(I)  with  polar  model  folded  in  rms  error  as  a  function  of 
sunspot  number  for  five  selected  polar  paths. 


6 


1  4 
RMS  ERROR 


2 


0 


Figure  18.  MINIMUF-85  with  polar  model  folded  in  rms  error  as  a  function  of 
sunspot  number  for  five  selected  polar  paths. 
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Although  the  size  of  the  polar  data  sample  was  small,  the  results  were 
unambiguous.  At  least  over  the  north  polar  paths  tested,  MINIMUF-3.5(I)  makes 
an  average  rms  error  in  MUF  prediction  of  10.92  MHz.  While  in  the  polar 
enhanced  versions  of  MINIMUF,  the  RMS  error  is  reduced  to  about  4.0  MHz.  It 
is  clear  that  a  considerable  improvement  in  predicting  MUF's  is  attained  by 
including  a  polar  model  in  the  MINIMUF  prediction  program. 
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DISCUSSION 


This  report  has  presented  the  development  of  MINIMUF-85,  a  simple  semi- 
empirical  algorithm  for  predicting  the  MUF.  It  was  developed  to  predict 
accurate  maximum  useable  frequencies  (MUFs)  under  conditions  of  anomalously 
high  sunspot  numbers,  to  predict  f qF2  values  suitable  for  ray-tracing  applica¬ 
tions,  and  to  predict  M  factor  values  useable  for  determining  the  mirror 
height  of  reflection  for  oblique  incidence  propagation,  and  to  improve  its 
accuracy  for  paths  into  or  crossing  the  polar  regions.  It  includes  sunspot 
number  dependence  in  both  the  fQF2  and  the  M  factor  calculations.  The  M 
factor  portion  of  the  algorithm  also  is  modified  to  include  seasonal  and 
diurnal  variations. 

This  section  of  the  report  discusses  several  topics  important  in  the 
application  of  MINIMUF-85  including  the  relationship  between  the  sunspot 
number  and  10.7-cm  solar  flux,  the  choice  of  solar  indices  for  forecasting, 
the  determination  of  take-off  angle  from  the  M-factor,  the  role  of  MINIMUF  in 
sounder  updating  of  schemes  for  determining  a  pseudo  sunspot  number,  and 
finally,  possible  future  improvements. 

RELATIONSHIP  BETWEEN  SUNSPOT  NUMBER 
AND  10.7-CM  SOLAR  FLUX 

It  has  been  common  practice  to  use  10.7-cm  solar  flux  as  a  measure  of 
solar  variation  in  MINIMUF  applications.  This  has  been  accomplished  using  the 
relationship  developed  by  Steward  and  Left in  (reference  26)  between  the  Ottawa 
10.7-cm  solar  radio  noise  flux  and  the  Zurich  smoothed  sunspot  number.  This 
relationship  is  given  by 

-4  2  (29) 

*12  *  63.7  +  0.728R12  +  8.9  x  10  \n 

where  4>  is  the  12-month  running  mean,  Ottawa  10.7-cm  solar  radio  noise  flux 
and  Rj2  is  the  12-month  running  mean  Zurich  sunspot  number.  It  was  determined 
using  data  from  November  1947  to  November  1968.  They  found  no  systematic 
differences  in  the  relationship  for  the  rising  and  declining  phases  of  the 
solar  cycles  investigated. 
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For  the  years  1947-1966,  Joachim  (reference  27)  presents  a  relationship 

between  R,_  and  *  ,  the  monthly  mean  of  the  daily  values  of  the  10.7-cm  solar 
li  m 

radio  flux  at  Ottawa.  It  is  given  by  the  expression 

-0.05R., 

♦  *  R10  +  46  +  23e  .  (30) 

m  12 

At  values  of  R12  between  10  and  80,  this  expression  gives  lower  values  for 
than  equation  (29)  does  for 

The  definition  of  R^  is 


fn+5 


[n-5 


in  which  Rk  is  the  mean  of  the  daily  sunspot  numbers  for  a  single  month  k  and 
R12  is  the  smoothed  index  for  the  month  represented  by  k  *  n. 

For  the  years  1947-1963,  Gladden  (reference  28)  presents  a  relationship 

between  Rz  ,  the  monthly  mean  Zurich  relative  sunspot  number,  and  #,  the 
m 

monthly  mean  daily  value  of  the  10.7-cm  solar  radio  flux  at  Ottawa.  Using  all 
of  the  data,  the  form  of  the  relation  was: 

*  =  0.9l(Rz  +  58.8  (32) 

m  ^  m) 

which  had  a  correlation  coefficient  of  0.98  and  a  rms  deviation  of  11.4. 

Since  the  solar  radio  emissions  in  the  range  of  wavelengths  from  3-cm  to 
30-cm  show  similar  variations  with  a  so-called  27-day  cycle,  Nicolet 
(reference  29)  determined  the  relationship  between  the  27-day  average  value  of 
the  Zurich  relative  sunspot  number  and  the  solar  radio  flux  at  10-cm  from  1951 
until  1962.  The  results  show  that,  for  R2?  >  50,  there  is  a  linear  relation¬ 
ship  between  these  quantities  which,  with  an  error  of  +10%,  takes  into  account 
the  variations  of  the  radio  flux  for  values  greater  than  $27  =  100  units, 
i.e. , 


$27  =  50  +  0.967  R27  . 

But  for  $27  <100  or  R27  <50,  the  linear  variation  is  given  by 


(33) 
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*27  -  68.0  +  0.607  R27  . 


(34) 


Thus,  the  close  association  of  the  10-cm  radiation  with  sunspots  cannot  be 
expressed  during  a  whole  cycle  with  the  same  linear  relationship.  Figure  19 
compares  the  variation  of  the  monthly  mean  and  27 -day  mean  $  provided  by 
equations  32,  33,  and  34  with  the  monthly  mean  and  27-day  mean  Zurich  sunspot 
number.  The  data  points  in  the  figure  are  for  the  period  February  1947  - 
December  1962.  The  nonlinear  curve  is  that  provided  by  equations  33  and  34. 
Note  that  the  relationship  provided  by  equations  33  and  34  appears  to  give  a 
better  fit  than  that  given  by  equation  32.  At  sunspot  zero  equations  30,  32, 
and  34  should  give  about  the  same  result  for  the  solar  flux  as  there  should  be 
very  little  difference  among  R^,  Rz^,  and  R^.  ^act»  only  equations  30 
and  34  give  comparable  results.  Hence,  equations  33  and  34  appear  to  provide 
a  more  accurate  estimate  of  the  mean  monthly  solar  10.7-cm  flux. 


VARIATION  OF  OTTAWA  10.7cm  FLUX  (<£)  WITH  ZURICH 
SUNSPOT  NUMBER  (RZ) 

(MONTHLY  MEANS. FEB.  1947-  OEC.ISSZ) 


Figure  19.  Variation  of  monthly  mean  and  27-day  mean  10.7-cm 
flux  with  monthly  mean  and  27-day  mean  Zurich  sunspot  number. 


CHOICE  OF  SOLAR  INDEX  FOR  FORECASTING 

MINIMUF-85,  like  MINIMUF-3.5  before  it,  will  be  used  to  estimate  a  near 
real-time  value  of  MUF.  Immediately  there  arises  :he  question  of  which  solar 
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index  to  use.  In  MINIMUF-85  two  indices  are  used.  In  the  polar  model  R^» 

the  smoothed  Zurich  sunspot  number  is  used.  In  the  rest  of  the  model  Rzm»  the 

monthly  mean  value  is  used.  Evaluation  of  equation  32  for  ♦  using  R^  an<i 

equations  33  and  34  for  <J>  using  Rzm  results  in  an  insignificant  difference 

in  *  for  either  R  going  from  0  to  300.  The  reason  for  this  is  that  the 

dispersion  of  the  monthly  mean  solar  flux  t  is  large  enough  to  cause  either 

expression  to  be  valid  (i.e.,  if  either  expression  is  used  to  estimate  the 

sunspot  number  from  $,  either  will  provide  equally  valid  results).  However, 

if  values  for  both  R,  „  and  Rz  are  available,  they  should  be  used  in  their 

12  m 

respective  parts  of  MINIMUF-85. 

The  dispersion  of  the  daily  values  of  10.7-cm  solar  flux  is  larger  than 
the  dispersion  in  the  27 -day  mean  values.  Figure  20  from  reference  29  shows  a 
plot  of  the  daily  values  for  1958  with  a  variation  reaching  ±20  percent  from 
equation  33.  The  figure  indicates  how  such  a  solar  index  may  lead  to 
important  errors  in  deducing  any  empirical  correlation  with  ionospheric 
parameters.  Gladden  (reference  28),  in  the  determination  of  a  relationship 
between  the  daily  values  of  the  solar  flux  *  and  Rz,  found  that  the  rms 

deviation  in  the  resulting  expression  for  was  32.6,  whereas,  the  rms 

deviation  for  this  expression  for  the  monthly  mean  was  only  11.4.  The 
resulting  expression  given  was 

$.  =  0.84  Rz.  +65.1  .  (35) 

d  d 

As  PROPHET  uses  equation  29  to  determine  the  sunspot  number  from  the  daily 
10.7-cm  flux,  the  results  of  using  these  two  expressions  is  compared  in  table 
21.  Note  that  the  use  of  equation  29  with  daily  flux  values  always  produces 
lower  sunspot  numbers  than  does  equation  35. 

To  further  examine  the  choice  of  solar  indices  to  be  used  with  MINIMUF- 
85,  we  compered  the  accuracy  of  using  the  various  choices  with  measured 
oblique  incidence  sounder  data  taken  during  Solid  Shield  exercises  conducted 
by  the  Naval  Research  Laboratory  which  took  place  between  3  May  and  20  May 

1981,  on  the  eastern  seaboard  of  the  United  States  (reference  31).  The 

locations  of  the  sounded  paths  are  indicated  in  figure  21. 
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Figure  20.  Dispersion  of  daily  values  of  1 0.7-cm  solar  flux  about 
27-day  mean  values  for  1 958. 


Table  21.  Comparison  of  two  expressions  relating  sunspot  number  R 
and  10.7-cm  daily  solar  flux. 


*12 

*d 

R 

63.7 

65.1 

0 

100.1 

107.1 

50 

145.4 

149.0 

100 

192.9 

191.1 

150 

244.9 

233.1 

200 

301.3 

275.0 

250 

362.2 

317.0 

300 
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Figure  21 .  Geometry  of  the  oblique  sounder  circuits  employed  by  NRL  during  Solid 
Shield  plotted  on  a  great  circle  map. 

Table  22  and  figures  22-25  show  the  nature  of  the  solar  and  geophysical 
conditions  during  the  Solid  Shield  exercises.  Table  22  shows  the  solar 
activity  background  for  each  of  the  days  during  May  1981.  It  gives  the  inter¬ 
national  sunspot  number  (the  continuation  of  the  Zurich  sunspot  number),  the 
3-  and  5-day  running  averages  of  this  number,  the  sunspot  number  Rz^  derived 
from  the  daily  10.7-cm  solar  flux  using  equation  35,  the  sunspot  number  R^ 
derived  from  the  daily  10.7-cm  solar  flux  using  equation  29,  the  daily  10.7-cm 
solar  flux  4^,  and  the  3-  and  5-day  running  averages  of  10.7-cm  solar  flux. 
Figure  22  shows  the  daily  variation  of  the  five  sunspot  numbers  given  in  table 
22.  Note  first  the  larger  variation  in  Rj,  the  international  sunspot  number 
than  in  either  the  running  daily  averages  and  those  produced  by  the  10.7-cm 
solar  flux.  Second,  note  that  values  produced  from  the  10.7-cm  solar  flux  are 
higher,  although  within  the  published  deviations  of  the  relationships,  than 
Rj.  Figure  23  shows  the  daily  10.7-cm  solar  flux  variation  given  in  table  22. 
Note  that  the  3-day  running  average  appears  between  the  other  curves  and  only 
lags  the  daily  curve  by  3  days.  Figure  24  is  a  plot  of  solar  flare  activity 
for  the  month  of  May  1981.  The  height  of  the  individual  bars  in  figure  24 
represents  the  relative  height  of  the  peak  of  each  solar  flare  in  energy  out- 
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Table  22.  Solar  activity  background  for  Solid  Shield,  May  1981. 


Day  SUNSPOT  NUMBERS  10.7-cm  SOLAR  FLUX 

of  3  and  5  3  and  5 


Month 

RI 

day  averages 

Rz . 
d 

R12 

*d 

day  averages 

Comment 

1 

112 

103 

91 

143 

142 

185.0 

177.3 

179.9 

2 

133 

112 

103 

149 

147 

190.4 

182.7 

179.6 

3 

156 

134 

120 

164 

160 

203.0 

192.8 

185.4 

burst 

4 

152 

147 

129 

181 

174 

217.5 

203.6 

193.7 

5 

162 

157 

143 

200 

189 

233.3 

217.9 

205.8 

6 

178 

164 

156 

193 

183 

227.1 

226.0 

214.3 

7 

171 

170 

164 

196 

186 

229.8 

230.0 

222.1 

8 

177 

175 

168 

182 

175 

218.3 

225.1 

225.2 

burst 

9 

158 

169 

169 

178 

171 

214.5 

220.9 

224.6 

10 

148 

161 

166 

177 

170 

213.4 

215.4 

220.6 

burst 

11 

169 

158 

165 

189 

180 
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Figure  25.  Magnetic  activity  for  May  1981  as  represented  by  the  index  Ap. 
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put  during  the  period.  The  solar  activity  was  high  during  the  full  term  of 
the  experimental  period  between  3  and  19  May  with  the  period  of  highest 
activity  between  7  and  14  May  1981.  During  the  month  of  May  1981  the  earth's 
magnetic  field  was  also  very  active.  This  is  shown  in  figure  25.  This 
activity  is  indicated  by  the  index  A^.  Note  that  until  9  May  1981  the 
magnetic  activity  was  very  low,  but  after  that  date  there  was  a  considerable 
increase. 

To  limit  the  amount  of  comparison  against  data  to  an  amount  handable 
(i.e.,  the  number  of  comparisons  are  proportional  to  the  number  of  indices 
being  considered)  data  for  5  and  7  May  1981,  Hurlbert-Norfolk  path,  were  used 
to  determine  the  accuracy  of  the  indices.  These  two  days  were  chosen  because 
they  occur  during  a  period  of  low  geomagnetic  activity.  Figures  26-29  show 
the  predicted  MUF  as  a  function  of  solar  index  for  each  of  the  two  days. 
Figures  26  and  27  show  the  effects  of  the  five  sunspot  indices  for  5  May  1981 
and  7  May  1981,  respectively.  On  both  days  the  sunspot  number  Rz  produced 
from  equation  35  using  daily  10.7-cm  solar  flux  4^  gave  the  highest  MUF;  the 
5-day  running  average  sunspot  number  gave  the  lowest  values.  Figures  28  and 
29  show  the  effect  of  using  different  10.7-cm  solar  flux  values  for  5  May 
1981,  and  7  May  1981,  respectively.  In  each  case  the  largest  values  were 
produced  by  the  dai'y  10.7-cm  solar  flux  4^.  Table  23  gives  the  bias  and  rms 
error  for  each  solar  index  for  the  two  days  during  the  Solid  Shield  exercise 
on  the  Hurlbert-Norfolk  path.  A  positive  bias  indicates  that  the  model 
predicts  high.  The  data  in  the  table  indicate  that  the  sunspot  number 
produced  from  the  daily  solar  flux  4^  using  equation  35  has  the  lowest  bias 
and  lowest  rms  error  for  the  sample  examined.  The  results  need  to  be 
confirmed  using  data  from  other  paths  and  other  exercises. 
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EFFECT  OF  SOLAR  ACTIVITY  FOR  SOLID  SHIELD 
HURLBERT-NORFOLK  MAY  5  188) 


5  DAY  SSN  AV. 


Figure  26.  Effect  of  sunspot  number  index  for  Solid  Shield,  Hurlbert-Norfolk,  5  May  1981 . 
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Figure  27.  Effect  of  sunspot  number  index  for  Solid  Shield,  Hurlbert-Norfolk,  7  May  1981 . 
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Figure  28.  Effect  of  10.7-cm  solar  flux  index  for  Solid  Shield,  Hurlbert-Norfolk,  5  May  1978. 
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Figure  29.  Effect  of  10.7-cm  solar  flux  index  for  Solid  Shield,  Hurlbert-Norfolk,  7  May  1981 
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Table  23.  Bias  and  rms  error  for  each  solar  index, 
Solid  Shield,  Hurlbert-Norfolk. 


Index  Value 

5  May  1981 

Bias  rms  error 

7  May  1981 

Bias  rms  error 

RI 

-t./L 

J  .  J7 

3-day  average 

-2.78 

3.63 

-2.84 

3.71 

5 -day  average 

-3.01 

3.78 

-2.92 

3.76 

Rzd 

-2.32 

3.37 

-2.58 

3.55 

R12 

-2.40 

3.42 

-2.67 

3.60 

3-day  10.7  cm  solar  flux 

-2.55 

3.49 

-2.67 

3.60 

5-day  10.7  cm  solar  flux 

-2.70 

3.58 

-2.76 

3.65 

SOUNDER  UPDATING 

The  day-to-day  variation  of  the  MOF  and  its  corresponding  predicted  MUF 
about  the  monthly  median  MOF  and  MUF,  respectively,  are  often  considerable  and 
may  exceed  the  combined  seasonal  and  solar  cycle  variation.  The  monthly 
median  of  MUF  reflects  the  main  effects  of  an  undisturbed  ionosphere.  The 
daily-hourly  MUF  deviates  from  the  median  by  2-  to  8-MHz  on  quiet  days  and  by 
much  larger  values  during  magnetic  storms.  The  causes  for  these  deviations 
are  complicated  and  cannot  be  compensated  for  by  simple  rules  of  thumb. 

The  purpose  of  sounder  updating  is  to  use  an  ionospheric  oblique 
incidence  sounder  to  provide  an  input  that  could  be  used  to  compensate  for  the 
more  systematic  part  of  the  deviation  of  MUF,  particularly  the  daily-hourly 
variation.  The  approach  used  is  to  obtain  a  maximum  observed  frequency  (MOF) 
on  a  reference  path,  and  use  this  information  to  determine  an  effective  or 
pseudo- sunspot  number  (references  31,  32,  33,  34,  35,  and  36).  This  pseudo¬ 
sunspot  number  is  then  used  in  MINIMUF  rather  than  a  daily  solar  index  to 
predict  the  MUF  on  other  paths  in  the  same  region. 

Figure  30  is  a  drawing  which  illustrates  the  approach  the  Naval  Research 
Laboratory  employs  to  perform  model  update  in  support  of  frequency  management. 
The  top  illustration  in  the  figure  indicates  the  nature  of  the  diurnal  varia¬ 
tion  of  both  the  measured  MOF  and  the  predicted  MUF  over  a  given  link.  There 
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is  a  bias  in  the  prediction  due  to  day-to-day  variability.  The  update  process 
is  shown  in  the  middle  portion  of  the  figure.  A  measurement  of  the  MOF  is 
obtained  from  an  oblique  sounding  over  a  known  circuit  and  input  into  the  up¬ 
date  algorithm.  The  algorithm  then  finds  the  sunspot  number  that  would  have 
produced  the  input  MOF  for  the  specific  time  of  the  measurement  over  the  path. 
This  is  accomplished  by  varying  the  sunspot  number  from  an  initial  value  of 
zero  until  the  condition  that  MOF  =  MUF  is  satisfied.  The  success  of  the  fit 
is  tested  by  determining  the  rms  error  over  a  24-hour  period  between  the 
predicted  MUFs  using  the  pseudo- sunspot  number  and  the  measured  MOFs  on  the 
same  path.  In  tactical  scenarios,  a  pseudo-sunspot  number  which  was  derived 
from  the  update  procedure  is  then  used  to  compute  MUFs  on  other  unknown  paths. 

The  model  update  is  deemed  successful  if  the  rms  error  of  the  experi¬ 
mental  paths  are  significantly  lower  than  that  yielded  by  employing  the  non- 
updated  model.  In  practice  it  has  been  found  that  updating  intervals  of  three 
hours  is  sufficient  and  that  rms  errors  of  1  MHz  in  some  cases  will  be 
obtained.  However,  values  of  2.5  MHz  are  more  common. 

To  compare  the  sounder  updating  technique  with  that  of  the  non-updated 
MINIMUF-85  results,  the  same  two  days  of  Solid  Shield  were  run  using  a  pseudo¬ 
sunspot  number  generated  for  2200  UT  (local  noon)  for  5  and  7  May  1981  as  used 
in  the  previous  section.  The  results  are  presented  in  table  24.  In  both 
cases  the  pseudo-sunspot  number  was  250.  A  slight  improvement  was  obtained 
over  that  of  non-updating.  The  sounder  updating  technique  is  limited  at  these 
high  sunspot  numbers  because  of  the  sunspot  number  limitation  built  into 
MINIMUF-85. 


Table  24.  Bias  and  rms  error  using  sounder  updated  sunspot 
numbers,  Solid  Shield,  Hurlbert-Norfolk. 


5  May  1981 

7  May  1981 

pseudo-sunspot  number 

”7515 

bias 

-2.16  MHz 

-2.39  MHz 

rms  error 

3.28  MHz 

3.43  MHz 
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EFFECTS  OF  UNDERLYING  LAYERS  ON  M- FACTOR  ESTIMATION 


In  the  calculation  of  MUF  in  MINIMUF-85,  the  MUF  is  determined  by  the 
product  of  the  fQF2  and  the  so-called  M-f actor.  The  range  dependence  in  the 
M-factor  is  based  on  the  work  of  Davies  (reference  20)  who  derived  a  family 
curves  for  the  M-factor  as  a  function  of  range  and  height  of  maximum  electron 
density.  The  range  dependence  used  was  obtained  by  curve-fitting  to  exact 
results  given  by  Davies  for  a  parabolic  layer  of  290  km  height  and  a  ratio  of 
semithickness/base  height  of  0.4.  This  corresponds  to  a  semi  thickness/height 
of  maximum  electron  density  of  0.29.  The  results  apply  to  a  single  parabolic 
layer  with  no  underlying  layers. 

Lockwood  has  presented  a  noniterative  procedure  which  enables  evaluation 

of  the  MUF  using  a  distance  factor  M(D)F2,  with  allowance  for  variations  in 

both  the  peak  height  and  changes  in  the  underlying  plasma  (reference  37).  The 

algorithm  presented  uses  values  of  the  ionospheric  characteristics  fQF2,  fQE 

and  M(3000)F2.  Use  is  made  of  the  Bradley-Dudeney  model  of  the  electron 

density  profile  (reference  38),  consisting  of  a  combination  of  linear  and 

parabolic  segments  to  represent  the  E-,  F1-,  and  F2-regions.  Lockwood  assumed 

a  fixed  value  for  the  ratio  Ym  F2/hm  F2  of  .29.  He  found  that  using  the 

M(3000)F2-  factor  scaled  according  to  standard  URSI  slider  resulted  in 

differences  between  values  obtained  by  using  ray-tracing  calculations  by  up  to 

7.5 7.  for  F2  peak  heights  less  than  500  km;  this  maximum  error  falling  to  3% 

for  h  F2,  the  F2  peak  height,  below  350  km.  However,  he  developed  a  correc- 
m 

tion  to  the  M( 3000 )F2-f actor,  scaled  from  the  ionogram,  based  on  the  ratio  of 

f  F2  to  f  E,  x,  that  was  always  accurate  to  within  0.5%.  Using  the  corrected 
o  o 

value  for  the  M(3000)F2-factor  in  the  algorithm  for  M(D)F2,  values  are 

obtained  with  accuracies  to  about  6%  for  any  range.  The  largest  errors  are 

for  ranges  in  excess  of  A000  km;  for  ranges  less  than  4000  km,  the  algorithm 

is  accurate  to  within  4%.  The  algorithm  presented  was  not  designed  for  use  at 

x  (f  F2/f  E)  less  than  1.95. 
o  o 

The  value  of  M(D)  can  be  calculated  in  terms  of  x,  M(3000)q,  and  D  by 

using 
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M(D) 


1  + 


[M(3000)q 


U 


(36) 


where  and  C^qqq  are  given  by  equations  37  and  38  below  for  the  ranges  D  and 
3000  km,  respectively.  The  parameter  is  given  by 

CD  -  0.72  -  0.628z  -  0.45z2  -  0.03z3  +  0.194z4  +  0.158z5  +  0.037z6  (37) 

where 


z  *  1 


2D 


D 

max 


The  value  of  the  maximum  range  in  kilometers  is  given  by 


D 

max 


3940  +  s 


1 

M(3000) 

o 


0.258 


(38) 


(39) 


with 

S  -  9900  +  mm 

X  X 


(40) 


Equations  36-40  allow  the  calculation  of  the  M-f actor.  This  is  done  in 
three  stages:  an  estimate  of  M(3000)q  at  a  range  of  3000  km  is  determined  in 
terms  of  M(3000)^  from  ionogram  data,  and  then  is  determined  at  a  range  of 
3000  km  and  finally  at  range  D  is  determined.  The  value  of  M(3000)q  is 
given  by 


M(3000)  =  M(3000).  -  0.124 

o  i 


+  [M(3000).z  -  4] 


0.0215 


+  O.OOSsin  [^^~L  -  1.9635^  . 


(41) 


The  value  of  M(3000)^  is  obtained  from  an  ionogram  or  equivalent  than 

obtained  from  a  numerical  map  representation  of  M(3000)^.  Alternatively,  it 

can  be  obtained  from  the  expression  given  by  Bradley  and  Dudeney  (reference 

38)  for  h  F2  if  h  F2  is  available: 
m  m 


h  F2 
m 


1490 


M(3000)F2  +  AM 


-  176 


(42) 
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where 


0.18  .  ,  , 

“  *7^X4  •  *  >  *-7  • 

For  no  underlying  layers  (i.e.,  x  *  100)  and  hmF2  *  290  km  (the  height 
used  currently  in  MINIMUF-85),  the  equations  become 

AM  =  1.82  x  10'3 

M(3000).  =  3.19560 

M(3000)  *  M(3000).  -  0.02 

o  1 

M(  3000 )q  =  3.17562 
s  -  9901.54 

D  =  4503.39  .  (A3) 

max 

The  value  of  M(D)  is  found  using  equations  36-38. 


The  assumption  given  by  Lockwood  for  the  ratio  hmF2/YmF2  of  3.5  implies 
for  the  Bradley-Dudency  profile  used  by  Lockwood  that  the  parameter  hmF2  can 
be  given  by 


h'F,  F2  +  104 


h  F2  = 
m 


f 0.613  )°’86 
[x-1.33j 


0.71429  + 


f  0.613  ^ 

{x  -  1.33J 


0.86] 


(44) 


where  h'F,F2  is  the  minimum  observed  virtual  height  of  reflection  from  the 
F- layer  or  F2-  layer  whichever  may  be  the  case.  When  x  =  100,  there  is  no 
underlying  ionization  and  hmF2  becomes 

h  F2  =  1.37564[h'F,F2  +  1.31598]  .  (45) 

m 

The  value  for  the  ratio  h  F2/Y  F2  of  3.5  also  implies  that  the  base  of  ioniza- 

m  m 

tion  of  the  layer  is  given  by 


h  F2  =  0.71429  h  F2  . 


(46) 


0  m 

If  there  were  no  ionization  below  the  F2-  layer,  h'F,F2  would  be  approximately 
equal  hQF2,  and  hmF2  would  be  approximately  equal  to  1.4  h'F,F2. 
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Contour  maps  or  numerical  coefficient  representations  of  h'F,F2  can  be 

used  to  infer  how  h  F2  differs  from  the  value  of  290  km  used  in  MINIMUF-85  or 
m 

even  be  used  to  determine  M(3000)F2  using  equations  42  and  44.  Figure  31, 
based  on  contour  maps  of  h'F,F2  (reference  39)  shows  how  h'F,F2  varies  as  a 
function  of  time  of  day,  season,  and  latitude.  The  F2-layer  reflection 
heights  given  in  reference  39  are  the  average  of  observations  for  the  years 
1944-1947.  Reference  39  also  shows  the  variation  of  F2-layer  height  with 
sunspot  number  for  June  at  1200  LT  at  four  receiver  sites.  There  is  generally 
an  increase  in  layer  height  with  increasing  sunspot  number  for  northern 
latitude  sites  and  a  decrease  with  increasing  sunspot  number  at  the  one 
southern  latitude  site.  Leftin  et  al.  (reference  40)  give  numerical  maps  of 
h'F,F2  base  height  as  a  function  of  season,  sunspot  number,  location,  and  time 
of  day. 
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Figure  31 .  Approximate  minimum  virtual  heights  of  F,F2-  layers,  January,  June,  December. 
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DETERMINATION  OF  TAKE-OFF  ANGLES 


In  the  evaluation  of  antenna  gains,  path  loss,  and  group  path  delay,  the 
elevation  is  necessary  to  determine  the  take-off  angles  of  the  possible  modes 
of  propagation.  This  is  done  by  calculating  the  so-called  mirror-reflection 
height.  This  is  the  height  at  which  an  equivalent  plane  mirror  would  have  to 
be  placed  to  reflect  unrefracted  waves  with  the  same  elevation  angles  at  the 
receiver  and  transmitted  as  for  the  real  ionosphere. 

In  PROPHET  the  take-off  angle  is  determined  using  figure  A. 10  (reference 
Al).  This  figure  illustrates  the  relationship  between  great  circle  distances 
in  km  and  radiation  or  elevation  angles  in  degrees  for  a  variety  of  propaga¬ 
tion  modes  at  fixed  heights. 

In  several  of  the  larger  HF  propagation  prediction  procedures,  the 
mirror-reflection  height  is  calculated  iteratively  (references  5  and  7); 
alternatively,  equations  based  on  some  mean  reference  model  ionosphere  are 
employed  (reference  A2).  The  advantages  of  an  iterative  process  are  that  the 
effects  of  variations  in  the  underlying  ionospheric  regions  and  of  spatial 
changes  can  be  taken  into  account  and  that  the  evaluation  can  be  continued 
until  the  required  accuracy  can  be  achieved.  For  many  applications,  however, 
it  is  undesirable  to  have  an  indeterminate  calculation  time,  and  the  accuracy 
with  which  the  ionization  density  profile  is  known  may  not  justify  such  a 
procedure. 

Lockwood  (reference  A3)  has  developed  new  algorithms  for  the  prediction 

of  the  mirror  reflection  height  hp.  Lockwood  (reference  AA)  has  shown  that 

the  evaluation  times  are  small  yet  the  accuracy  is  higher  than  that  of  many 

existing  simplified  procedures.  The  algorithms  allow  for  the  effects  of 

variations  in  underlying  ionization  but  neglect  the  effects  of  the  geomagnetic 

field  and  apply  for  a  horizontally  stratified  ionosphere.  A  fixed  ratio  of 

3.5  for  h  F2/y  F2  was  assumed  in  his  development.  Lockwood  used  the  procedure 
m  m 

for  determining  the  M(D)F2  factor  described  in  reference  37  to  calculate  ttp  at 
a  range  D  and  fixed  r  (f/MUF)  by  iterating  the  elevation  angle.  In  this  way 
hp(D)  curves  are  calculated  for  various  r  and  for  each  of  the  model  profiles. 
Over  a  large  part  of  the  range  of  r,  hp  at  a  fixed  D  is  only  weakly  dependent 
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on  r  (for  r  ■  0.8  to  0.95  when  h  F2  ■  500  km,  x  *  3.33  and  for  r  =  0.7  to  0.95 

m 

when  h  F2  ■  250  km,  x  =  3.33,  except  at  the  shortest  ranges).  In  every  case 
m 

where  r  >  0.7,  tVj,  is  greatest  at  all  D  when  r  is  unity;  it  is  much  greater 
than  for  the  lower  r  at  the  smallest  and  largest  distances.  Based  on  these 
variation  patterns,  Lockwood  developed  two  representations  of  the  mirror- 
reflection  height:  1)  at  r  ■  1.0  for  propagation  at  the  MUF  or  above;  2)  at  r 
between  0.75  and  0.95  for  propagation  at  the  frequency  of  optimum  transmitting 
(FOT).  Because  no  formulas  are  presented  for  lower  ratios  of  r,  the  latter  is 
often  used  at  lower  values  of  the  ratio. 


At  all  but  the  shortest  ranges,  hp  can  be  described  by  the  linear  form 
when  0.75  <  r  <  0.95 

h_  =  sD  +  C  .  (47) 

For  M(3000).  <3.57 


s 


1 

M(3000)i 


0.24 


B  . 


(48) 


For  M(3000)i  >  3.57 

s  »  0.04B  .  (49) 

The  parameter  B  varies  linearly  with  the  inverse  of  the  fourth  power  of  x: 

B  «  0.03  +  14.0/x4  .  (50) 

The  parameter  C  is  given  by 

c  -  358  +  {roW;  -  °-35]  <51) 

where  a  is  given  by 

a  =  1880  -  3200/x5  .  (52) 

Equations  47-52  can  be  combined  into  the  equation 

Kj,  =  358  -  (11-100  a)^18.8  -  +  aD^0.03  +  |^Jkm  (53) 
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where  a  ■  [1 /Ml  3000^  -  0.24]  or  0.04,  whichever  is  the  larger. 

Fo.  no  underlying  ionization  (x  *  100)  and  a  hmF2  of  290  km,  becomes 

Kj,  -  290.61  +  0.0022D  km  .  (5/») 

The  general  form  of  the  h^D)  curves  for  r  =  1.0  cannot  be  approximated 
by  a  single  linear  relationship  as  is  possible  for  the  0.75  <  r  <  0.95  curves. 
In  this  case  the  fit  is  given  by 

IVj,  —  SjW  +  c^  +  4  w  <  0.95  (55) 

Kj.  =  thjlw  =  0.95  w  >  0.95 

where  w  equals  D/D  and 
max 


The  intercept  c^  is  given  by 


-  35  + 


r  i 

al [m(3 


M(3000) 


-  0.225 


where 

Oj  *  1785  -  4000/x3  . 

The  slope  s^  varies  linearly  with  M(3000)q  to  the  power  -1.5: 

Sl  =  230  +  B.  M(3000)  _1*5  -  0.14 
1  1  o 


(58) 


(59) 


=  325  +  6.4  x  10A/x3'8  (60) 

Equations  55-60  can  be  used  to  evaluate  Hp  for  f  =  MUF  at  a  given  D,  in 
coniunction  with  values  for  M(3000)  and  obtained  in  the  evaluation  of 

J  O  IUclX 

the  M(D)  factor  from  the  previous  section. 
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For  no  underlying  layer  and  an  hmF2  of  290  km,  these  equations  become 
hj,  -  195.5  +  0.0537  D  +  23^A503>4  -  lj  .  (61) 

At  D  *  1000  km,  hj,  «  329.8  km  for  f  -  MUF;  and  hj,  *  292.8  km  for  0.75  <  r  < 
0.95. 


For  realistic  values  of  h^,,  w  used  in  equations  55  and  56  should  have  a 
minimum  value  given  by 

500-h  F2 

w  ,  =  o.i  -  0.025  ■  ~Trr~~  •  (62) 

min  310 

The  problem  with  the  use  of  equation  53  for  calculating  hj,  for  0.75  <  r  <  0.95 

at  short  range  is  that  h^,  is  no  longer  linear  with  decreasing  range.  At  short 

range  the  h^,  variation  is  similar  to  that  for  r  ■  1.0  (frequencies  near  the 

MUF),  only  not  as  pronounced.  That  is  as  r  decreases  from  1.0  to  0.6,  the 

3 

increase  in  hj,  with  decreasing  distance  D  is  proportional  to  r  .  Further,  as 
x  decreases  from  a  value  of  10.0  to  2.0,  the  increase  in  h^  at  short  ranges 
becomes  less  significant.  A  correction  factor  is  found  to  be  applied  to 
equation  53  for  ranges  less  than  This  correction  is  found  by  finding 

the  range  for  which  equation  55  had  a  minimum.  This  minimum  range  is  inserted 
in  equation  55  for  r  *  1  to  determine  h^  Then  Ah^,  is  found  from  Ah^  = 

•v  •  Yhin-  dmin  is  *iv“  ‘y 


D  =  D  /  — 
MIN  MAX  Sj 


(63) 


and  Afyj,  is  given  by 

AhT  =  23  DMAx(d  '  D^)  +  D^{D  '  Dmin)  *  (64) 

The  restriction  of  equation  62  still  applies.  The  final  correction  A'Kj, 
is  given  by 
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(65) 


The  angle  of  transmission  from  earth  (sometimes  called  takeoff  angle), 
written  in  terms  of  the  mirror  reflection  height  hj,  and  the  great-circle 
distance  D  from  transmitter  to  receiver  is 


(66) 


where  n  =  number  of  hops  and  re  =  radius  of  the  earth. 


FUTURE  IMPROVEMENTS  IN  MINIMUF 

Before  termination  of  the  discussion  section  of  this  report,  possible 
additional  improvements  that  might  be  made  to  MINIMUF  are  discussed.  There 
are  three  areas  where  improvements  might  be  made.  These  are  in  tve  fQF2 
model,  in  the  M- factor  model,  and  in  the  Polar  region  model. 


fQF2  Representation 


The  representation  of  fQF2  in  MINIMUF- 85  is  given  by 

f  F22  =  A  +  A.  cos  1/2  xoff 
o  o  I  ell 

where  is  a  function  of  sunspot  number  and  represents  the  amplitude  varia¬ 
tion  of  f  F2  (at  x  «  =  0°.  A  +  A.  =  maximum  diurnal  value  of  f  F2  and  at 

0  6X1  o  i  2  ^ 

Y  =  90°.  A  =  minimum  diurnal  value  of  f  F2  ).  If  examined  the  contour 
Aeff  o  o 

plots  of  f  F2  given  in  references  21  and  22  show  that  the  minimum  and  maximum 

values  of  fQF2  depend  on  sunspot  number,  geographic  location,  and  season.  By 

using  the  f  F2  data  base  used  to  get  the  representation  for  A.  now  in  MINIMUF- 
o  ^  1 

85,  it  should  be  possible  to  obtain  Aq,  A1  and  cos  'xgff  as  a  function  of  sun¬ 
spot  number,  season,  and  geographic  location.  At  each  fQF2  measurement  site 

A  ,  A,  ,  and  n  would  be  determined  for  each  season  and  sunspot  number  that 
o  1 

minimizes  the  diurnal  error  in  f^F2.  Having  determined  A^,  A^,  and  r)  at  each 
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site,  then  the  functional  form  of  geographic  location,  sunspot  number,  and 
season  dependences  of  these  parameters  can  be  determined. 

M- factor  Representation 

The  M-factor  representation  could  be  improved  by  introducing  the  effects 
of  the  underlying  layers  on  M-factor  estimation  given  earlier  by  equations 
36-41.  The  parameter  M(3000)i  can  be  found  from  equation  42  using  equation  44 
for  hmF2.  The  parameter  h'F,F2  required  by  equation  44  can  be  determined 
either  from  numerical  maps  of  h'F,F2  (reference  39)  or  from  figure  31.  To 
fine  tune  the  model,  the  resultant  MINIMUF  model  will  need  te  be  fit  to  an 
oblique  sounder  data  much  as  the  present  version  was. 

Polar  Representation 

Possible  improvements  in  the  polar  model  include 

•  development  of  a  polar  M-factor  model, 

•  introduction  of  a  seasonal  dependence  into  the  folding  function  to 
reduce  the  winter  bias  error,  and 

•  refinement  of  the  adjustable  parameters  to  further  minimize  errors. 
Finally,  south  Polar  MOF  data  is  needed  to  test  the  model  in  that 
region. 
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CONCLUSIONS 


We  have  described  several  changes  to  the  MUF  algorithm  used  in  the 
PROPHET  family  of  prediction  programs.  The  main  effect  of  the  changes  is  to 
improve  the  capability  of  the  model  during  periods  of  very  high  solar 
activity,  a  situation  where  the  previous  model  had  known  deficiencies. 

We  have  also  improved  the  fQF2  model  to  provide  a  much  more  accurate 
prediction  capability  for  this  parameter.  This  will  improve  short  distance 
MUF  predictions  and  give  more  accurate  fQF2  values  for  use  in  other  PROPHET 
outputs,  such  as  ray  tracing,  where  this  parameter  is  required  as  input. 

We  have  changed  the  M-factor  model  to  include  sunspot  numbers,  seasonal, 
and  diurnal  dependencies.  These  changes  allow  more  flexibility  in  using  this 
parameter  for  other  PROPHET  output,  such  as  a  prediction  ly  which  requires 
the  parameter  as  an  input. 

Lastly,  we  have  introduced  a  special  fQF2  for  use  in  the  polar  regions  of 
the  world.  This  is  accomplished  by  the  use  of  a  folding  function  which  folded 
in  the  polar  f  F2  model  and  folded  out  the  MINIMUF-85  fQF2  model.  The  folding 
function  is  a  function  of  geomagnetic  latitude  and  sunspot  numbers.  Inclusion 
of  the  polar  model  in  MINIMUF-85  reduced  the  overall  bias  in  the  model  from 
0.16  MHz  low  to  0.14  MHz  low  and  reduced  the  rms  error  from  4.19  MHz  to  4.08 
MHz.  On  five  Polar  paths  it  changed  the  bias  and  rms  erroi,  respectively, 
from  6.85  MHz  low  and  10.92  MHz  to  C.6  MHz  low  and  4.0  MHz. 
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APPENDIX  A 


BASIC  PROGRAM  FOR  MINIM1TF-85 

The  listing  of  MINIMUF-85  that  follows  is  written  in  extended  BASIC  for 
the  Tektronix  4052A  computer.  The  subroutine  Razgc  determines  the  latitude 
and  longitude  of  a  point  on  a  great  circle  path  given  the  range  and  bearing 
from  the  point.  The  subroutine  Gcraz  gives  the  range  and  bearing  between  two 
points. 


INPUT  PARAMETERS 


Z3  =  transmitter  latitude. 


radians 


(f  <  Z3  <  f) 


Z4  =  transmitter  longitude,  radians  ( - 2tt  <  Z4  <  2tt) 
Z5  *  receiver  latitude,  radians  <  Z5  < 


Z6  =  receiver  longitude,  radians  (-2n  <  Z6  <  2it) 

Z0(1)  =  path  length,  radians 

Z0(2)  *  azimuth  of  path,  radians 

Ml  =  month 

DO  =  day  of  month 

HO  =  hours,  GMT 

MO  =  minutes,  GMT 

S9  =  monthly  median  sunspot  number 

PI  =  3.141593 

PO  =  1.5707963268 
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OUTPUT  PARAMETER 


The  output  is  as  follows: 
J9  -  MUF  in  MHz. 
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100  !  MNIMUF85  TEST  "MNIMUFTEST" 

110  ! 

120  !  INPUT  DATA  DRIVER  FOR  MINIMUF 

130  Rl-PI/180 

140  P1-2*PI 

150  PO-PI/2 

160  DIM  ZO (8) 

170  Z3-75*R1  !  TRANSMITTER  LATITUDE 

180  Z4*-125*R1  !  TRANSMITTER  LONGITUDE  WEST  POSITIVE 

190  Z5=51.95*R1  !  RECEIVER  LATITUDE 

200  Z6= 1 76 . 58*R1  !  RECEIVER  LONGITUDE  WEST  POSITIVE 

210  Ml=l  !  MONTH  1-  JAN 

220  D0= 1 5  !  DAY 

230  !H0«0  HOUR 

240  M0=0  !  MINUTE 

250  S9=9  !  SMOOTH  SUNSPOT  NUMBER 

260  A*Z5 

270  B=Z6 

280  C-Z3 

290  D-Z4 

300  CALL  Gcraz 

310  ZO ( 1 ) =R  !  RANGE  BETWEEN  TRANSMITTER-RECEIVER 
320  Z0(2)-S  !  BEARING  ANGLE 
330  FOR  H0=0  TO  23  STEP  1  !  LOOP  FOR  24  HOURS 

340  CALL  Mnimuf 

350  PRINT  " LAT=" ; Z3/R1 , " LON= " ; Z4/R1 , " LAT»" ; Z5/R1 , " LON= " ; Z6/R1 

360  PRINT  " MONTH= " ;M1 , "DAY=" ;D0, "HOUR=" ;H0, "MINUTE=" ; MO 

370  PRINT  "SSN*";S9 

380  PRINT 

390  PRINT  "MUF=" ;J9 

400  NEXT  HO 

410  END 

420  ! 

430  SUB  Mnimuf 

440  DIM  Ssn(6) ,Scn{6) 

450  Saa»0.814*S9+22.23 

460  Sab=l. 3022-0. 001 56*S9 

470  FOR  J-l  TO  6 
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480  Sar-2*J*PI*M1/12 

490  Ssn(J)-SIN(Sar) 

500  Scn(J)-C0$(Sar) 

510  NEXT  J 

520  Sac-0 . 9925+0 . 01 l*Ssn ( 1 ) +0 . 087*Scn ( 1 )  -0 . 043*Ssn ( 2 ) +0 . 003*Scn ( 2 ) 

530  Sac-Sac-0. 013*Ssn(3)-0.022*Scn(3)+0.003*Ssn(4)+0.005*Ssn  (5) 

540  Sac-Sac+0.018*Scn(6) 

550  I-H0+M0/60 

560  J-INT (Z0( 1 )/0 . 62784)+! 

570  L-1/(2*J) 

580  Ak6-1.59*Z0(l) 

590  IF  Ak6<l  THEN 

600  Ak6-1 

610  END  IF 

620  J9-100 

630  IF  Z0(1)>0. 94174  THEN 

640  Ak-2*J-1 

650  ELSE 

660  Ak-J 

670  END  IF 

680  Kkk-1/Ak6 

690  IF  Kkkol  THEN 

700  Kkk-0.5 

710  END  IF 

720  FOR  D9-1  TO  Ak  STEP  1 

730  IF  Z0(1)>0. 94174  THEN 

740  Akl-D9*L 

750  ELSE 

760  Akl»l/(2*Ak6)+(D9-l)*(0.9999-l/Ak6) 

770  END  IF 

780  C-Akl*Z0(l) 

790  D-Z0(2) 

800  A-Z5 

810  B-Z6 

820  CALL  Razgc 

830  E-R 

840  F-S 

850  IF  F«»0  THEN  870 

860  F-F+Pl 

870  G-0 . 0172*( 10+(M1 - 1 )*30 . 4+DO) 
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880  Slt«I-F/(Rl*15) 

890  IF  SI t<24  THEN  910 

900  Slt-Slt-24 

910  IF  SI t>0  THEN  930 

920  Slt-Slt+24 

930  A-E 

940  B-F 

950  Smg«0.9792*SIN(A)+0.2028*C0S(A)*C0S(B- 1.2043) 

960  Sdg-ASN(Smg) 

970  IF  ABS(Sdg)<0. 95993  THEN 

980  Sgf-0 

990  ELSE 

1 000  Sgf-0 . 7578*SQR  ( l+3*Smg*Smg )  *0 . 5  -  0 . 5 

1010  END  IF 

1020  H«0.409*C0S(G) 

1030  P«3.82*F+12+0.13*(SIN(G)+1.2*SIN(2*G}) 

1040  IF  P<«24  THEN  1070 

1050  P-P-24 

1060  GO  TO  1090 

1070  IF  P«>0  THEN  1090 

1080  P-P+24 

1090  Q-2.5*Z0(l)*Kkk  MIN  PO 

1100  Q-SIN(Q) 

1110  Q-1+2.5*Q*SQR(Q) 

1120  IF  COS (E+H) >-0 . 26  THEN  1170 

1130  G-0 

1140  S-0 

1150  Sad-1 

1160  GO  TO  1610 

1170  S«(-0.26+SIN(H)*SIN(E))/(C0S(H)*C0S(E)+l .0E-3) 

1180  S-S  MAX  -1  MIN  1 

1190  S-12-ASN(S)*7.6394 

1200  T-P-S/2 

1210  IF  T«>0  THEN  1230 

1220  T-T+24 

1230  U-P+S/2 

1240  IF  U<=24  THEN  1260 

1250  IMJ-24 

1260  V*ABS(COS(E+H) ) 

1270  W-9.7*VA9.6 
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1280  W-W  MAX  0.1 

1290  X-I 

1300  IF  U<T  AND  (I-U)*(T-I)>0  THEN  1440 

1310  IF  U»>T  AND  ( I-T)*(U- I )<-0  THEN  1440 

1320  IF  T<-I  THEN  1340 

1330  X-X+24 

1340  Y-PI*(X-T)/S 

1350  Z-PI*W/S 

1360  F-(T-X)/W 

1370  F-F  MAX  -100  MIN  100 

1380  G«V*(SIN(Y)+Z*(EXP(F)-COS(Y)))/(l+ZA2) 

1390  P-V*(Z*(EXP(-S/W)+l))*EXP((S-24)/2)/(l+ZA2) 

1400  IF  G->P  THEN  1420 

1410  G-P 

1420  Sad-l.ll-0.01*Slt 

1430  GO  TO  1610 

1440  IF  U<-I  THEN  1460 

1450  X-X+24 

1460  Stt-X-U 

1470  Stu-14*Stt/(24-S) 

1480  Sag-PI*(Stu+l)/15 

1490  Sah-2*Sag 

1500  Sal-1. 0195-0. 06*SIN(Sah)-0.037*COS(Sah)+0.018*SIN(2*Sah) 

1510  Saj--0.003*COS(2*Sah)+0.025*SIN(3*Sah)+0.018*COS(3*Sah) 

1520  Sak-0 . 007*S IN (4*Sah ) - 0 . 005*C0S  (4*Sah )  +0 . 006*S I N ( 5*Sah ) 

1530  Sal -0 . 017*C0S(5*Sah) -0 . 009*SIN(6*Sah )  -0 . 004*COS (6*Sah ) 

1540  Sad-Sal+Saj+Sak+Sal 

1550  Z-PI*W/S 

1560  F-(U-X)/2 

1570  F-F  MAX  -100  MIN  100 

1580  Y-S/W 

1590  Y-Y  MAX  -100  MIN  100 

1600  G-V*(Z*(EXP(Y)+1))*EXP(F)/(1+ZA2) 

1610  Sae-Sab*Sac*Sad 

1620  H«SQR(6+Saa*SQR(G) )+Sgf 

1630  H»H*(l-0.1*EXP((S-24)/3)) 

1640  H-H*(1+(1-SGN(Z5)*SGN(Z3))*0.1) 

1650  H-H*(l-0. 1*{1+SGN(ABS(SIN(A) ) -COS(A) ) ) ) 

1660  CALL  Fof2 

1670  J9-J9  MIN  H 
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168C  PRINT  "H  -  ",H 

1690  NEXT  D9 

1700  J9-J9  MAX  2  MIN  50 

1710  END  SUB 
1720  SUB  Fof2 

1730  ! LOCAL  Gg,Ff,Ys,Za,Am,Phi ,Tmo,Rltm,Rlgm,Plr,Bb,T,U,V,Y,Z,C,W,X 

1740  Phi-ST t*PI/12 

1750  Tmo-Ml+(D0+I/24)/30-0.5 

1760  Rltm-Sdg 

1770  Rlgm«C0S(A)*SIN(B-1.2043)/C0S(Rltm) 

1780  Rlgm-Rlgm  MAX  -1  MIN  1 

1790  Rlgm-ASN(Rlgm) 

1800  X=(2.2+(0.2+S9/1000)*SIN(Rltm))*COS(Rltm) 

1810  Ff-EXP(-(XA6) ) 

1820  Gg=l-Ff 

1830  T»PI*Tmo/12 

1840  V-SIN(T) 

1850  U«COS(T+T) 

1860  Y=SIN(Rlgm/2) 

1870  Ys*COS(Rl gm/2- PI/20) 

1880  Z-SIN(Rlgm) 

1890  Za-SQR(ABS(Z) ) 

1900  Am-l+V 

1910  IF  Sdg<0  THEN  1960 

1920  C*-23.5*PI/180 

1930  W=EXP(- 1 . 2*(COS(R1  tm+C*C0S(Phi ) )  -COS  (R1  tm) ) ) 

1940  Plr«(2+1 .2*S9/100)*W*(1+Q.3*V) 

1950  GO  TO  2000 

1960  B=V*(0.5*Y-0.5*Z-YA8)-Am*U*(Z/Za)*EXP(-4*Y*Y) 

1970  PI r«2 . 5+2*S9/100+U*(0 . 5+(  1 . 3+0 . 2*S9/100)*YsA4) 

1980  PI r*Pl r+( 1 . 3+0 . 5*S9/100)*COS( Phi  - PI*( 1+B) ) 

1990  Plr=Plr*(l+0.4*(l-V*V))*EXP(-l*V*YsA4) 

2000  T =Gg*HA2/8 .12+0. 66*F  f *P1  r 

2010  IF  T<0  THEN  2040 

2020  Ff2=TA0.5*2.85 

2030  H-Ff2*Q*Sae 

2040  END  SUB 

2050  SUB  Gcraz 

2060  IF  ABS (A-C) >1  .OE-5  OR  ABS ( B-D) >1 . 0E-5  THEN  2100 
2070  R-1.0E-6 


2080  S-0 

2090  GO  TO  2290 

2100  IF  AB$(A-P0)>l.OE-5  THEN  2140 

2110  R-PO-C 

2120  S-PI 

2130  GO  TO  2290 

2140  IF  ABS(A+P0)>1.0E-5  THEN  2180 

2150  R-PO+C 

2160  S-0 

2170  GO  TO  2290 

2180  E-SIN(A) 

2190  F-COS(A) 

2200  G-SIN(C) 

2210  H-E*G+F*COS(C)*CO$(B-D) 

2220  H-H  MAX  -1  MIN  1 

2230  R-ACS(H) 

2240  S«(G-E*H)/(F*SIN(R)) 

2250  S-S  MAX  -1  MIN  1 

2260  S-ACS(S) 

2270  IF  SIN(B-D) >0  THEN  2290 

2280  S-Pl-S 

2290  END  SUB 
2300  SUB  Razgc 

2310  IF  Ol.OE-5  THEN  2350 

2320  R-A 

2330  S-B 

2340  GO  TO  2610 

2350  IF  ABS(A- PO) >1  .OE-5  THEN  2390 

2360  R-PO-C 

2370  S-B 

2380  GO  TO  2610 

2390  IF  ABS (A+PO) >1 . OE-5  THEN  2430 

2400  R-C-PO 

2410  S-B 

2420  GO  TO  2610 

2430  E-SIN(A) 

2440  F-COS(A) 

2450  G-COS(C) 

2460  H-E*G+F*SIN(C)*COS(D) 

2470  H-H  MAX  -1  MIN  1 
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2480  P-ACS(H) 

2490  Q«(G-E*H)/(F*SIN(P)) 

2500  Q=Q  MAX  -1  MIN  1 

2510  Q-ACS(Q) 

2520  R-PO-P 

2530  IF  S IN (D) >0  THEN  2560 

2540  S-B+Q 

2550  GO  TO  2570 

2560  S-B-Q 

2570  IF  SoPl  THEN  2590 

2580  S-S-Pl 

2590  IF  S=>-P1  THEN  2610 

2600  S=S+P1 

2610  END  SUB 
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APPENDIX  B 


FORTRAN  PROGRAM  FOR  MINIMUF-85 

The  listing  of  MINIMUF-85  that  follows  is  written  in  FORTRAN  77  for  the 
HP  9050  computer.  The  parameters  passed  to  the  subroutine  and  those  returned 
by  it  are  described  in  the  comments  portion  of  the  routine. 
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cp 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

Cl 


c 

cm 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

Cl 


subroutine  muf85  (tlat,  tlon,  rlat,  rlon,  itime,  cpnt ,  9in,cmuf ) 
subroutine  muf85 
update  nov  1985 

an  improved  verilon  of  muf35  which  includes  ssn, season  and 
diurnal  dependence  in  the  M-factor  plus  an  improved  F0F2  model 
call  muf85(tlat  ,tlon,rlat,rlon, it i me, ssn , cmuf ) 

this  routine  computes  the  maximum  useable  frequency  (cmuf)  for 
a  given  propagation  path.  the  required  input  is: 

parameters  passed. 

tlat  -  transmitter  latitude  in  radians 
tlon  -  transmitter  west  lingitude  in  radians 
rlat  -  receiver  latitude  in  radians 
rlon  -  receiver  west  longitude  in  radians 
itime  -  six  element  array  containing  the  month,  day 
hour,  minute,  Julian  day,  and  year 
ssn  -  sunspot  number 

parameters  returned: 

cpnt  -  path  control  point  info  in  radians 
cmuf  -  classical  muf  in  megahertz. 

called  by  subroutine  or  function:  mufluf 

subroutines  and  functions  called:  fof£ 

path 

razgc 

sygn 

common  blocks  referenced:  hite 

logical  v,  first 

int  eger  i  t i me( 6  > 

real  cpnt(8),  k5,  10,  lmt,  k8,  k9,  m9,  mlat , sn( 6 ) , cn( 6 ) 

common  /hits/  h,v,ym 
common  /hite/  h,  v,  ym 

h  is  the  height  of  the  path 

v  is  a  logical  variable  which  decides  if  the 
path  length  is  calculated  from  1000  km  f 
end  points  and  if  the  multilayered  ionospheric 
model  Is  to  be  used.  these  will  be  performed 
i  f  v  i s  t  rue . 

ym  is  the  f  layer  thickness. 


data  pi/3. Ml  59265/, tuopi/6 .2831 853/, half pi/1 .57079632/, 

4  dtr/0. 01  7453293/, rtd/57. 2957795/, r0/6371 ./ 

data  s8  /250.O/,  fml  /0.728/,  fms  /0. 52998/, 

4  f m2  /0. 00356/,  fm3  /63.75Z,  fm4  /0. 00178/, 

4  f i rst  / . t  rue . / 
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1 5  ■  float(  itim#(3)  )  +  floatl  itime(4)  )/60.0 
c 

c  convert  10.7  cm  flux  to  sunspot  numbsr 

c  ssn  ■  (  sqrt (  fms  -  fm2#(  fm3  -  fluxlO  )  )  ~  fml  )/fm4 

c  ssn  *  amax1(  am;in1(  ssn,  250.0  ),  0.0  ) 

c 

c  determine  number  of  hops 

c  1  hop  for  path  length  <•  4000  km  (  0.6278  radians  ) 
c  2  hop  otherwise 

c  a3  is  a  6th  order  fourier  series  based  on  month  which  is  part 

c  of  the  new  M-factor 

c  a2  is  a  linear  function  of  ssn  in  the  M-factor 

c  al  is  a  linear  function  of  ssn  in  the  critical  frquency 

c  expression 

c 

do  500  n  =  1,6 

qn  *  float(2*n) 

arg  *  pi *qn* i t im*( 1 )/1  2 . 0 

sn(n)  *  sin(arg) 

cn(n)  *  cos(arg) 

500  continue 

a3  *  . 9925+ . 01 1 *sn( 1 )+ . 087*cn( 1 )- . 043*sn(2) 

1  + . 003*cn ( 2  )  - . 013*sn(3)-.022*cn(3) 

2  +. 003*sn(4)+. 00S*sn(5>+. 018*cn(6) 
al  «  .  81  4*ssn+22 . 23 

a£  a  1 .3022- . 00156*ssn 
call  path( 1 1  at , t Ion , r 1  at , r Ion , cpnt ) 
gl  »  cpnt ( 1 ) 
azim  *  cpnt ( 8  ) 
c 

c  control  point  changes  for  muf85 

c 

i f ( . not .  v Ithen 
h3*1 .59 

ak6  *  h3*cpnt ( 1  ) 

if (ak6  . It  .  1  .  0 )akG«1  . 0 

k5  *  1  .  0 / a k 6 

if(k5  ne.  1.0)  k5  ■  .5 

khop  ■  int ( cpnt ( 1 )/. 62784 )+1 

kkhop  =  khbp 

i  f  ( cpnt  ( 1  )  .gt.  0. 94174)kkhop  *>  2*khop-1 

els* 

c 

c  old  control  point  methor  for  raytrac* 

c 

khop  =  1 

if  (  gl  .gt.  0.62784  )  khop  =  2 
k5  ■  1 . 0/float (  khop  ) 
kkhop  =  khop 
end  i  f 
c 

emuf «1 00,0 
ym  *  100.0 
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do  160  k 1  *  1 , k  khop 
c 

c  10,  wO  ■  latitude  and  west  longitude  of  control  points 
c  mid-point  for  1  hop  case;  points  2000  km  from  each 
c  end  for  2  hop  case  . 
c 

if  (  khop  .eq.  1  )  pi  “  gl/2.0 

if  (  khop  eq.  2  and.  kl  eq.  1  )  pi  °  0.31392 

if  (  khop  .eq.  2  .and.  kl  .eq.  2  )  pi  -  gl  -  0.31392 

c 

c  if  V  is  .false.,  do  cntrl  pt  calculations  like  in  apes 
c 

if(v)  go  to  600 
c 
c 

c  control  point  method  for  muf85 

c 

if(cpnt<1)  .gt.  .94174)  then 
x  k 1  =  kl 
xhop  *  khop 
akl  **  x  k  1  /  (  2 . 0*  xhop  ) 

else 

akl  *  1 . 0/(2 . 0*ak6)+float ( k1-1)*( . 9999-1 . 0/ak6) 
end  i  f 
c 

pi  *  gl *ak 1 

c 

600  call  razgc(  rlat,  rlon,  pi,  azim,  10,  wO  ) 
c 

c  lmt  *  local  mean  time  in  hours  at  the  control  point 
c  mlat  -  geomagnetic  latitude  at  the  control  point 
c 

if  (  wO  .  ge  0.0  )  then 
lmt  =  wO 

else 

lmt  =  wO  +  twopi 
end  i  f 

lmt  *  t £  -  lmt *rtd/1 5 . 0 
if  (  lmt  It.  0.0  )  t  hen 
lmt  *>  lmt  +  24.0 
else  if  (  lmt  .  ge  .  24.0  )  then 
lmt  3  lmt  -  24.0 
end  i  f 

smg  *  0.9792*sin(  10  )  +  0.2028*cos(  10  )*cos(  wO  -  1.2043  ) 
smg  *  amaxH  aminl  (  smg,  1.0  >,  -1.0  ) 
mlat  »  asin(  smg  ) 
c 

c  gyro  frequency  for  lat  >  55deg 

c 

if  (  abs(  mlat  )  It.  0.95993  )  then 
gyro  3  0.0 
else 

gyro  =  0 . 3789*sqrt (  1.0  +  3.0*smg*smg)  -  0.5 
end  i  f 

c 

c  yl  =  2*p i *date/365 . 25 
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c  y2  ■  -#olar  d#clination 
c  k8  ■  tim*  of  local  noon 
c 

y1»0. 0172*1  10.0  +  float (itimell 1-1  1*30.4  ♦  itim#<2)  ) 
y2»0 . 409*cos( y  1  ) 
c 

k8«*3 . 82*w0+ 1 2 . 0+0.13*(«in(y1 1  +  1 .2*#in(2.0*y11) 
if  1  k8  .gt.  24.0  )  than 
k 8  -  k8  -  24. 0 

•la#  if  1  k8  .1#.  0.0  1  th#n 
k8  -  k8  +  24.0 
•nd  i  f 
c 

c  m9  «  m-factor  =  muf/f0f2 
c 

m9  *  aainf (  2.5*g1*k5,  halfpi  1 
m9  *  ainl  m9  ) 
m 9  ■  1.0  +  2.5*m9*«qrt<  »9  > 

c 

c  Chang##  to  includ#  altitud#  of  th#  f-lay#r  variation# 
c  for  diurnal , lat itud#,  and  #olar  cycl# 
c  (if  v  is  .falsa.  bypa«*  all  this  variation  #tuff> 
c 

i f 1 . not . v )  go  to  50 

51  cchi»si n( 1 0 ) *s i n ( -y2 ) +cos 1 1 0 ) *co*( -y2 )*co#(  (t5-k8)*15.*dtr) 

chi*acoa(cchi ) 
c 

c  altitud#  variation#  of  f-lay#r 
c 

chi  1 »( ( chi  +  0 . 349 ) / . 873 ) **2 
xl-( 10/0.5241**2 

xmax«3 . - ( ##n-25 .  )*5.#-3+1 .25*co#( (  10- 0.961*0. 045) 
d*lx-2 . *co«(chi-0 .873+.  698*co#<  (  10-0.961*0 . 045)  1 
6  +( s#n-25 .  )*0.01*#xp(-chi1)*#xp(-xl) 

x11*xmax  -  d*lx 

c 

50  continue 

if  1  cost  10  +  y2  )  gt.  -0.26  1  go  to  100 
c 

c  no  daylight  on  path  at  any  tim#  during  th#  day 
c 

gO  ■  0.0 
k9  «  0.0 
go  to  140 
100  continu# 
c 

c  k9  ■  l#ngth  of  daylight 
c  t  ■  tim#  of  #unris# 

c  1 4  ■  tim#  of  #un##t 

c 

k9  *  (  -0.26  +  #in( y2 1 *#in(  1 0 1  )/(  co# ( y2 ) *co# ( 1 0 1  +  1  . 0#-3  1 
k9  =  amaxH  aminl  (  k9,  1.0  1,  -1.0  1 
k9  *  12.0  -  asinl  k9  1*7.6394 
t  -  k8  -  k9/2 . 0 

if  (t.  It.  0.0  1  t  ■  t  +24.0 
t4  »  k8  +  k9/2 . 0 
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if  (  t 4  .  gt  .  S4.0  )  t4  “  t4  -  24.0 
cO  ■  abs 1  cos (  10  +  y2  )  ) 
t9»9 . 7* 1 amax 1 1 cO ,  .  1  ) ) +*9 . 6 
1 9  •  amax 1  (  t9 ,  0.1  ) 

1 6  -  t5 

if  (  (  t4  .It.  t  .and.  1t5-t4)*(t-tS>  . gt .  0.0  )  or. 

!  <  1 4  .ge.  t  .and  ( t5-t ) * ( t4-tS )  .la.  0.0  >  )  go  to  120 

c 

c  day  time  at  control  point 
c 

if  (  t  . gt .  t5  )  t6  *  t6  +  24.0 
c 

c  local  time  conversion 

c  t5  is  local  time,  wO  is  longitude  in  radians 

c 

z  *  w0*(180 . 0/3. 14159265) 
c 

c  local  time  dependant  factor  for  M-factor 

c 

hrlcl  ■  t5  -  z/15.0 

i f (hrlcl  ge.  24.0)hrlcl  <*  hrlcl  -  24.0 
iflhrlcl  .It.  0.0)  hrlcl  »  hrlcl  +  24.0 
a4  =  1.11-.01  *  hrlcl 

c 

g9  «  pi  * (  t6  -  t  )/k9 
g8  *  p i *t  9/k  9 
u  *  (  t  -  t6  )/t9 

u  *  aminl (  amaxl (  u,  -87.0  >,  +87.0  ) 
u 1  »  -k9/t  9 

ul  ■  aminl  (  amaxlt  ul  ,  -87.0  ),  +87.0  ) 
g0»c0*(sin(g9)+g8+<exp(u>-cos(g9) ) )/( 1 . 0+g8*g8) 
gl  »  c0*(  g8*t  *xp(  ul  )  +  1.0  )  ) 

4  *exp (  (  k9  -  24.0  )/2.0  )/t  1.0  +  g8*g8  ) 

gO  *  amax 1 (  gO ,  g3  ) 

if(v)  yin  »  100.0*<1  .  0  +  0 . 2*cos ( pi  * ( 1 1 6-t ) /k9  -  0.5))) 
go  to  140 
1 £0  cont i nue 

c 

c  night  time  at  control  point 

c 

c 

c  6th  order  fourier  series  night  time  factor  for  M-factor, 

c  based  on  hours  after  sunset, t2 

c 

if<t4.gt.t5)  t6  *  t6  +  24.0 
1 1  ■  t  6  -  t  4 
t2  *  1 4 . 0  +  t 1 / 1 24 . 0-k 9 ) 
ag  ■  pi  * ( 1 2+ 1 . 0 ) /I  5 . 0 
agl  *  2.0*ag 
ag2  =  4.0*ag 
ag3  *  6 . 0*ag 
ag4  ■  8 . 0*ag 
ag5  ■  1 0 . 0*ag 
ag6  =  1 2 . 0*ag 

a14  -  1.0195  - . 06*sin<ag1 )-. 037*cos< agl )+ . 01 8*sin(ag2) 
a24  *  - . 003*cos1ag2 )+ . 025*sin(ag3)+ . 018*cost ag3) 
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*34  >  . 007*sin(ag4 )- .  005*cos(ag4)+ . 006*sin(ag5) 

*44  ■  . 01 7+cos(ag5 )- . 009*sin( ag6 )- . 004*cos(ag6 ) 

*4  ■  al 4+a24+a34+a44 

c 

g8“pi*t9/k9 
u  -  (  1 4-t 6 )/2 . 0 

u  *  aminl (  amax1(  u,  -75.0  ),  +75.0  ) 
u 1  -  -k  9/t  9 

u  1  *  aminl  (  amax1(  ul ,  -75.0  ),  +75.0  ) 
g0*cO* ( g8* ( explul  )+1.0))*#xp(u)/(1 . 0+g8*g8) 

140  continu* 

i f ( . not . v  )  go  to  150 

h  *  aminl (350.0, amax1<250.  ,x1  1  *ym>) 
c  th*  slop#  of  th#  mfactor  variation 

xm«-1  .*-3* am ini (6 .0*g1/(khop*.31  ) ,6.  ) 
c  the  new  mfactor 

m9*m9+xm* ( h-£90 ) 
c 

c  g2  ■  muf  at  control  point 
c 

150  continue 

g2  -  sqrt (6.0  +a1  *  sqrt(gO))  +  gyro 

g2  *  g£*(  1.0  -  0.1*#*p(  (  k9  -  24.0  )/3.0  )  ) 

g2  ■  g£*(  1.0  +  (  1.0  -  sygn(  tlat  )*sygn<  rlat  )  )*0.1  ) 

g2  *  g£* (  1.0  -  0 . 1*<  1.0  +  sygn(  abs(  #in(  10  )  ) 

*  -  cos (  10)  )  )  > 

if  (  abs(  mlat  )  .g*.  0.95993  )  th#n 

c 

c  F0F2  corrects  for  plar  region  F0F2.  result  is  C2  if 

c  not  in  polar  regiond 

c 
c 

g2  ■  m9*f of£(  g2,  lmt,  itimt,  10,  wO,  mlat,  ssn  ) 

else 

g£  *»  g2*m9 
end  if 
c 

g2  a  g£  *  a2*a3+a4 

cmuf  «  aminl (  cmuf ,  g£  ) 

160  continu* 

cmuf  ■  aminl  (  amax1(  cmuf,  2.0  ),  50.0  ) 
c 

return 

end 
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c 


function  fof£<  ff2.  lmt,  itime,  lat,  Ion,  ml«t,  ««n  ) 

function  fof2 

x  *  f of 2 (  f  f  2 ,  lmt,  itimt,  Ut ,  Ion,  mlat,  sen  ) 


reference  (to  ba  supplied  whan  available) 


input  ■ 
f  f2 
lmt 
it  ima 

lat 

Ion 

mlat 

isn 


critical  frequancy  from  muf35  in  mhi 
local  maan  tima  at  lat, Ion  in  hours 
integer  array  containing  month,  day, 
hour,  minute,  Julian  day,  and  year 
geographic  latitude  in  radians 
geographic  west  longitude  in  radians 
magnetic  latitude  in  radians 
sunspot  number 


- • real 

-  real 

-  integer 

-  real 

-  real 

-  real 

-  real 


output  '■ 
f  of  2 


the  f 2-layer  critical  frequancy  in  mhz  -  real 


called  by  subroutine  or  function:  muf35 


subroutines  and  functions  called:  none 


common  blocks  referenced: 


none 


integer  itime<6) 

real  lat,  lmt,  Ion,  mlat,  mlon 

data  Pi  /3 • 1415926/ 


c 


fli  “  mmeulVt  it  ime  (  2  )  ♦  it  ime(  3 )  /84  0  +  i  t  ime  C  4 )  / 1  440 . 0 

&  -0.5 

cmlat  ■  cos(  mlat  ) 

mlon  =  cos(  lat  )*sin<  ion  -  1  . 2043  l/cmi at 
mlon  =  amax  1  (  aminH  mlon,  1 
mlon  *  asin(  mlon  ) 
x  =  (  2.2  +  (  0.2  +  ssn/1000 

ff  ■  exp (  -(  x**&  >  5 


),  -10  ) 

)*sin(  mlat  )  )*cmlat 


gg  ■  1.0  -  f  f 
t  =  pi*tmo/1 2 . 
v  *  a  i  n  (  t  ) 
if  (  mlat  .ge. 
w  =  exp<  -1  • 
plr  =  (  2.0 


0 


0.0  )  then 

£*(  cos(  mlat  -  0.41015+cos(  phi 
+  0 . 0 1 2*ssn  )*w*(  10  +  0 . 3*v  ) 


cmlat  ) 


else 

u  B  cos (  tat  ) 

y  =  s i n (  ml  on/2 . 0  ) 

ys  »  cos(  m lon/2.0  -  pi/20  0  ) 

i  «  sin(  mlon  ) 

xs  »  sqrt (  abs (  z  )  ) 


)/30 . 0 


) 
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am  ■  1.0  +  v 

b  *  v*(  (  y  -  i  )/£.0  -  y**8  )  -  am*u*(  z/za  )*axp( 
ya4  ■  ys**4 

plr  »  £.5  +  ssn/50.0  +  u*(  0.5  +  (  1.3  +  0.00£*ssn 
6  +  (  1.3  +  0.005*ssn  )*coa(  phi  -  pi*(  1.0 

6  +  (  1.0  +  0.4*1  1.0  -  v*v  )  )*exp(  -v*ys4 

and  if 

f  of  £  ■  2 . 85*«qrt (  gg*f f£*f f£/8 .  1  £  +  0.66*ff*plr  ) 

raturn 

and 


-4.0*y*y  ) 

)*ys4  ) 

+  b  )  ) 

) 
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cp 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

cz 

c 

c 

c 

c 

c 

c 

c 


c 

r 

c 

c 

c 

c 


subroutine  path  (tlat,  tlon,  rlat,  rlon,  cpnt ) 
subroutine  path 

call  p'th(tlat,tlon, rlat, rlon, cpnt) 

this  routine  computes  the  range ,  az  iinut  h ,  and  control  point 
coordinates  for  a  given  propagation  path.  the  method  assumes 
a  spherical  earth  with  a  radius  of  6371-  km.  the  required 
input  for  this  module  is: 

tlat  transmitter  latitude  in  radians 

tlon  transmitter  west  longitude  in  radians 

rlat  receiver  latitude  in  radians 

rlon  receiver  west  longitude  in  radians 

this  subroutine  returns  the  following  information  in  an  8  word 
real  array  ( cpnt  )  : 

cpnt ( 1  )  distance  between  the  receiver  and  transmitter  in 
rad i ans 

cpnt(£)  latitude  of  midpoint  in  radians 
cnpt(3)  west  longitude  in  radians 

cpnt ( 4 )  latitude  of  point  1000km  from  the  receiver  in  radians 
cpnt(5)  west  longitude  of  point  1000km  from  receiver  in 
rad i ans 

c p n t ( 6 )  latitude  of  point  1000km  from  transmitter  in  radians 
cpnt(7)  west  longitude  of  point  1000km  from  transmitter 
in  radians 

cpnt(8)  azimuth  from  receiver  to  transmitter  in  radians 

cpnt(4)  through  cpnt(7)  will  not  be  computed  for  paths  less  than 
1000  km  (0.15696  radians)  in  length. 

subroutines  and  functions  used:  gcraz 

razgc 

common  blocks'  none 


d i mens i on  cpnt ( 8 ) 
get  range  and  azimuth 

call  gcraz(  rlat,  rlon,  tlat,  tlon,  cpnt ( 1  ) ,  cpnt(S)  ) 
get  mid-point  coordinates 
pi  =  cpnt ( 1  )  /£  .  0 

call  razgc(  rlat,  rlon,  pi,  cpnt(B),  cpnt(2),  cpnt(3)  ) 
is  path  length  >=  1000  km? 

if  (  cpnt ( 1 )  It.  0.15696  )  go  to  100 
yes  -  get  coordinates  of  1000  km  points 
pi  =  0. 15696 
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call  razgc(  rlat,  rlon,  pi,  cpnt(8),  cpnt(4)(>  cpnt(5)  ) 
pi  ■  cpnt(l)  -  0.15696 

call  razgc(  rlat,  rlon,  pi,  cpnt(8),  cpnt ( 6 ) ,  cpnt(7>  ) 
100  continue 
return 
end 
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cp 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

Cl 

c 


subroutine  razgcl  latl,  lonl,  range,  aziin,  latS,  long  > 
•ubroutine  razgc 

call  r azgc ( 1 at  1  ,  lonl , range  , az i m ,  lat2, 1  on2 ) 

this  routine  computes  the  latitude  and  west  longitude 
( i«t2,  long)  of  a  point  a  specified  range  from  a  given 
point  on  the  earth's  surface,  also  required  for  input 
is  the  azimuth  (azim)  to  the  new  point  in  radians.  this 
method  assumes  a  spherical  earth  and  recognizes  the 
degenerate  cases  of  the  given  point  being  at  the  north 
or  south  pole.  for  the  degenerate  cases,  azim  should  be  0 
or  pi  and  long  is  undefined.  however,  azim  is  not  checked, 
and  long  is  arbitrarily  set  equal  to  lonl.  this  routine 
recognizer,  the  degenerate  case  when  range  is  set  to  zero, 
all  coordinates  are  in  radians 

subroutines  and  functions  used'-  none 

common  blocks:  none 


real  latl,  lonl,  lat2,  long 


data  pi/3. 14159/,  twopi/6 .88318/, halfpi/t .570796/ 
data  rt d/57. 295779/, dtr/0. 01 7453/ 
c 

c  test  for  degenerate  cases 

if  (  abs<  latl  -  halfpi  )  .  gt  .  1 . Oe-5  )  go  to  100 


c 

c  the  given  point  is  the  north  pole 
c 

latS  3  halfpi  -  range 
long  *  lonl 
go  to  £ 0 0 
1 0  0  cont 1 nue 

if  (  abs(  latl  +  halfpi  )  gt  .  1.0e-5  )  go  to  120 


c 

c  the  given  point  i3  the  south  pole 


c 

lat2  =  range  -  halfpi 
long  *  lonl 
go  to  200 

120 

cont inue 

if  (  range  ,  gt  .  0.0  ) 

go  to  130 

c 

c 

point  £  coincident  with 

point  1 

c 

1  at 2  3  latl 
long  *  lonl 

go  to  800 
130  continue 
c 
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c  general  case 
c 

si  ■  sin<  latl  ) 
cl  *  cos(  latl  ) 
c2  ■  cos(  range  ) 

ca  -  s1*c2  +  c1*sin(  range  )*cos(  azim  ) 
ca  «  a*ln1(  amaxU  ca,  -1.0  ),  +1.0  ) 
a  =  acos(  ca  ) 
c 

c  test  if  destination  ends  up  on  the  poles 
c 

if(  abs ( a ) . gt . 1 . 0e-5  )  go  to  140 
lat£  »  halfpi 
lon£  *  lonl 
go  to  £00 
140  continue 

if(  abs(a-pi)  gt .  1 , 0e-5  )  go  to  150 

1  at 2  «  -halfpi 
lon2  ■  lonl 
go  to  £00 
1 50  cont i nue 
c 

c  everything  seems  ok,  get  destination  coordinates 
c 

eg  *  (  c£  -  si *ca  )/<  c1+sin(  a  )  ) 
eg  *  aminK  amaxH  eg,  -1.0  ),  +1.0  > 
g  ■  aco*(  eg  ) 
lat£  «  halfpi  -  a 
sa  ■  sin(  azim  ) 

if  (  sa  .ge.  0.0  )  long  *  amod(  lonl  -  g,  twopi  ) 

if  (  sa  It.  0.0  )  lon£  =*  amod(  lonl  +  g,  twopi  ) 

200  continue 
return 
end 
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function  sygn  (  y  ) 
cp 

c  real  function  sygn 

c 

c  x«sygn(y) 

c 

c  this  function  returns  the  value  of  0  if  y  is  0,  -1.  if  y  is 

c  less  than  zero  and  a  +1  .  if  y  is  greater  than  zero, 

c 

c  subroutines  and  functions  used:  none 

c 


c 

c 

cz 

c 

common  blocks : 

none 

if  (y)  100,  £00, 

300 

1  00 

sygn  *  -1.0 
go  to  999 

200 

sygn  *  0.0 
go  to  999 

300 

c 

sygn  =  1.0 

999 

return 

end 
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